Tutorial¶
This tutorial is also available as jupyter and org-mode notebook.
Quick start¶
Berni provides a database of interaction models and trajectory samples to quickly start molecular dynamics or Monte Carlo simulations. The focus is on simple models to study liquids and glasses. The package strives to follow FAIR-data practises.
The package provides an API to browse and retrieve models and samples:
models: database of all available modelsmodels.get("lennard_jones"): get a specific model by its qualified namesamples: database all available samplessamples.get("lennard_jones-13ce47602b259f7802e89e23ffd57f19.xyz"): get a copy of a sample by its qualified name
Both models and samples are slightly customized TinyDB databases, which you can query using the TinyDB methods.
The interaction models can be exported in formats suitable for simulation with a range of backends, such as atooms, RUMD, and LAMMPS.
Interaction models¶
Let’s see how to browse the models database.
Available models¶
The metadata of the available models can browsed by iterating over berni.models
import berni
for model in berni.models:
print(model['name'])
bernu_hiwatari_hansen_pastore
coslovich_pastore
dellavalle_gazzillo_frattini_pastore
gaussian_core
grigera_cavagna_giardina_parisi
harmonic_spheres
kob_andersen
kob_andersen_2
lennard_jones
roux_barrat_hansen
roy_heyde_heuer-II
roy_heyde_heuer
wahnstrom
These are the fields (or columns) of the database
berni.models.fields()
['_model', 'description', 'doi', 'name', 'notes', 'reference']
You can pretty print the database like this
berni.models.pprint(['name', 'description'])
name description
--------------------------------------------------------------------------------------------------------
bernu_hiwatari_hansen_pastore Binary soft-sphere mixture with size ratio of 1.4
coslovich_pastore Short-ranged pairwise-additive model for silica
dellavalle_gazzillo_frattini_pastore Binary Lennard-Jones mixture model for NiY alloys
gaussian_core One-component Gaussian core model with long-range cutoff
grigera_cavagna_giardina_parisi Binary soft-sphere mixture with size ratio of 1.2 and smooth cutoff
harmonic_spheres Binary mixture of harmonic spheres with size ratio of 1.4
kob_andersen Binary Kob-Andersen Lennard-Jones mixture
kob_andersen_2 Ternary Kob-Andersen Lennard-Jones mixture
lennard_jones One-component Lennard-Jones model
roux_barrat_hansen Binary soft-sphere mixture with size ratio of 1.2
roy_heyde_heuer-II variant of roy_heyde_heuer with optimized parameters
roy_heyde_heuer 2d model of a silica bilayer
wahnstrom Binary Lennard-Jones mixture with size ratio of 1.2
Since berni.models inherits from TinyDB, you can use all the TinyDB methods to browse and search the database, see the TinyDB docs.
The name field is a readable but “fully qualified” name that you can use to retrieve the specifications of the interaction model
Get the one-component Lennard-Jones model and inspect the parameters
model = berni.models.get('lennard_jones')
pprint(berni.models.get("lennard_jones"))
{'cutoff': [{'parameters': {'rcut': [[2.5]]}, 'type': 'cut_shift'}],
'potential': [{'parameters': {'epsilon': [[1.0]], 'sigma': [[1.0]]},
'type': 'lennard_jones'}]}
The model specifications are in a dict / json format, with two possibile schemas (1 and 2). The default schema requires the potential and cutoff fields. Additional fields are stored in the metadata field, but are not returned by the get method.
Here is how the default schema looks like
pprint(berni.schemas[1])
{'$schema': 'https://json-schema.org/draft/2020-12/schema',
'properties': {'cutoff': {'items': {'properties': {'parameters': {'type': 'object'},
'type': {'type': 'string'}},
'required': ['type', 'parameters'],
'type': 'object'},
'type': 'array'},
'metadata': {'properties': {'doi': {'type': 'string'},
'name': {'type': 'string'},
'notes': {'type': 'string'},
'reference': {'type': 'string'}},
'required': ['name'],
'type': 'object'},
'potential': {'items': {'properties': {'parameters': {'type': 'object'},
'type': {'type': 'string'}},
'required': ['type', 'parameters'],
'type': 'object'},
'type': 'array'}},
'required': ['potential', 'cutoff'],
'type': 'object'}
Schema version 2 stores the interaction parameters by pairs of species, instead of a single array.
pprint(berni.models.get("kob_andersen", schema_version=2))
{'potential': [{'cutoff': {'parameters': {'1-1': {'rcut': 2.5},
'1-2': {'rcut': 2.0},
'2-2': {'rcut': 2.2}},
'type': 'cut_shift'},
'parameters': {'1-1': {'epsilon': 1.0, 'sigma': 1.0},
'1-2': {'epsilon': 1.5, 'sigma': 0.8},
'2-2': {'epsilon': 0.5, 'sigma': 0.88}},
'type': 'lennard_jones'}]}
This is useful to generate interaction potentials for some third-party backends.
To choose or build your model, we need some more details about the schema layout, so keep reading.
Potentials and cutoffs¶
Each model is associated to a list of interaction potentials along with their parameters. If we limit ourselves to two-body potentials, we can write
where \(\alpha\) and \(\beta\) are chemical species in a mixture.
Turning back to the LJ model
model['potential']
[{'type': 'lennard_jones', 'parameters': {'epsilon': [[1.0]], 'sigma': [[1.0]]}}]
Notice how the epsilon and sigma parameters are both 2d arrays, as they depend on the chemical species of the interacting particles. There is a cutoff field as well, since the potentials are typically cut off at some distance.
model['cutoff']
[{'type': 'cut_shift', 'parameters': {'rcut': [[2.5]]}}]
The potentials and cutoffs have standardized names
for potential in sorted(berni.potentials):
print(potential)
fene
gaussian
harmonic
inverse_power
lennard_jones
sum_inverse_power
yukawa
The cutoffs have standardized names and some aliases too
for cutoff in sorted(berni.cutoffs):
print(cutoff)
cubic_spline
cut
cut_shift
linear
linear_cut_shift
quadratic
quadratic_cut_shift
shift
Trajectory samples¶
berni provides a database of trajectory samples to make it easy to
start a new simulation
check the implementation of interactions in third-party code
provide a little repository for further analysis
Each sample is associated to one the models defined above. Of course, you are free to use these samples for other interaction models.
Browsing the database¶
Browse the database of trajectory samples.
berni.samples.pprint(['name', 'model'], sort_by='model')
name model
----------------------------------------------------------------------------------------------------
bernu_hiwatari_hansen_pastore-f61d7e58b9656cf9640f6e5754441930.xyz bernu_hiwatari_hansen_pastore
coslovich_pastore-488db481cdac35e599922a26129c3e35.xyz coslovich_pastore
grigera_cavagna_giardina_parisi-0ac97fa8c69c320e48bd1fca80855e8a.xyz grigera_cavagna_giardina_parisi
kob_andersen-8f4a9fe755e5c1966c10b50c9a53e6bf.xyz kob_andersen
roy_heyde_heuer-II-b8d70742799933357ea83314590d2b4d.xyz roy_heyde_heuer-II
lennard_jones-5cc3b80bc415fa5c262e83410ca65779.xyz lennard_jones
Get a local copy of a Lennard-Jones fluid sample.
local_file = berni.samples.get("lennard_jones-13ce47602b259f7802e89e23ffd57f19.xyz")
By default, local_file will be stored in a temporary directory. You can use it to start a simulation, benchmarking your code or further analysis using other packages.
We rely on TinyDB to perform queries on the metadata of the available samples. For instance, here we look for samples of the Lennard-Jones models with a unit density
query = berni.Query()
samples = berni.samples.search((query.model == 'lennard_jones') &
(query.density == 1.0))
pprint(samples)
[{'density': 1.0,
'format': 'xyz',
'md5_hash': '13ce47602b259f7802e89e23ffd57f19',
'model': 'lennard_jones',
'name': 'lennard_jones-13ce47602b259f7802e89e23ffd57f19.xyz',
'number_of_particles': 256,
'potential_energy_per_particle': -3.8079776291909284,
'url': 'https://framagit.org/coslo/berni/-/raw/master/berni/samples/lennard_jones-13ce47602b259f7802e89e23ffd57f19.xyz'},
{'density': 1.0,
'format': 'xyz',
'md5_hash': '5cc3b80bc415fa5c262e83410ca65779',
'model': 'lennard_jones',
'name': 'lennard_jones-5cc3b80bc415fa5c262e83410ca65779.xyz',
'number_of_particles': 108,
'potential_energy_per_particle': 0.0,
'temperature': 0.0}]
We can tidy up the output with pprint()
berni.samples.pprint(['name', 'density', 'number_of_particles'],
cond=(query.model == 'lennard_jones') & (query.density == 1.0))
name density number_of_particles
------------------------------------------------------------------------------
lennard_jones-13ce47602b259f7802e89e23ffd57f19.xyz 1.0 256
lennard_jones-5cc3b80bc415fa5c262e83410ca65779.xyz 1.0 108
Here is another way to inspect the full metadata of a given sample, identified this time by the file name
query = berni.Query()
pprint(berni.samples.search(query.name == "lennard_jones-5cc3b80bc415fa5c262e83410ca65779.xyz"))
[{'density': 1.0,
'format': 'xyz',
'md5_hash': '5cc3b80bc415fa5c262e83410ca65779',
'model': 'lennard_jones',
'name': 'lennard_jones-5cc3b80bc415fa5c262e83410ca65779.xyz',
'number_of_particles': 108,
'potential_energy_per_particle': 0.0,
'temperature': 0.0}]
For more advanced queries, check out what TinyDB can do!
The database may provide some provenance information via the notes optional field. We do not aim at “full reproducibility” here.
Exporting models for backends¶
In addition to specifying the interaction models, berni can generate potentials objects for different interaction “backends”. At present there there is support for atooms, LAMMPS and RUMD - more to come!
Atooms¶
The model specifications of berni are compatible with the native f90 atooms backend. Here we compute the potential energy of the configuration stored in local_file using lennard_jones model.
from atooms.backends.f90.trajectory import Trajectory
from atooms.backends.f90 import Interaction
with Trajectory(local_file) as th:
system = th[0]
model = berni.get("lennard_jones")
system.interaction = Interaction(model)
print(system.potential_energy(per_particle=True))
-3.807977629191845
We can compare it to the reference value stored in the samples database and check they agree. This is useful if you want to test a new implementation of a backend.
query = berni.Query()
sample = berni.samples.search(query.name == "lennard_jones-13ce47602b259f7802e89e23ffd57f19.xyz")[0]
sample["potential_energy_per_particle"]
-3.8079776291909284
RUMD¶
In addition to the native f90 backend, there is currently support for RUMD. It is possible to generate potential objects for RUMD like this
potentials = berni.export("kob_andersen", "rumd")
The potentials list can then be passed to an Simulation instance defined in RUMD to carry out a simulation. You may want to use the atooms.rumd wrapper for that.
LAMMPS¶
Berni will try to export the model using the “analytical” potential functions defined in LAMMPS and fallback to tabulated potentials in all the other cases
import berni
berni.export("kob_andersen", 'lammps')
pair_style lj/cut 1.0
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2
pair_modify shift yes
berni.export("bernu_hiwatari_hansen_pastore", "lammps")
pair_style table spline 10000
pair_coeff 1 1 /tmp/tmpast7x_o3/phi.1-1 POTENTIAL
pair_coeff 1 2 /tmp/tmpast7x_o3/phi.1-2 POTENTIAL
pair_coeff 2 2 /tmp/tmpast7x_o3/phi.2-2 POTENTIAL
The tabulated potential files are generated automatically.