API Reference

dompap

Simulations of point-like particles in any dimension with any pair potential. The package uses NumPy and Numba for efficient calculations and SymPy to implement any pair potentials. The user is not expected to be familiar with these packages but only basic Python syntax.

Repository: <urpedersen/dompap>

class dompap.Simulation(particle_types: ~numpy.ndarray = <factory>, positions: ~numpy.ndarray = <factory>, box_vectors: ~numpy.ndarray = <factory>, image_positions: ~numpy.ndarray = <factory>, velocities: ~numpy.ndarray = <factory>, betas: ~numpy.ndarray = <factory>, neighbor_list: ~numpy.ndarray = None, positions_neighbour_list: ~numpy.ndarray = <factory>, masses: ~numpy.ndarray = <factory>, pair_potential_r_cut: ~numpy.float64 = np.float64(1.0), neighbor_list_skin: ~numpy.float64 = np.float64(1.0), max_number_of_neighbors: ~numpy.float64 = np.int32(512), number_of_neighbor_list_updates: int = 0, number_of_steps: int = 0, time_step: ~numpy.float64 = np.float64(0.01), temperature_target: ~numpy.float64 = np.float64(1.0), temperature_damping_time: ~numpy.float64 = np.float64(0.1), pair_potential: callable = None, pair_potential_str: str = None, pair_force: callable = None, pair_curvature: callable = None, epsilon_func: callable = <factory>, sigma_func: callable = <factory>)

The Simulation class

The default simulation is a 5x5x5 simple cubic lattice of particles with unit mass, diameter and epsilon. The box vectors are [5, 5, 5] and the particles are placed at [0, 0, 0], [0, 0, 1], … [4, 4, 4]. The default pair potential is the harmonic repulsive potential (1-r)**2, where r is the distance between two particles. The default pair potential is truncated at r_cut = 1.0. The default neighbor list is updated if a particle has moved more than the skin distance, which is 0.5 by default. The default maximum number of neighbors is 128. The ntegrator is a Langevin VT Leap-frog integrator. The default time step is 0.01, the default target temperature is 1.0, and the default temperature damping time is 0.1. Below is an example of how to set up a simulatio (using the default values):

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_positions(unit_cell_coordinates=([0.0, 0.0, 0.0],), cells=(5, 5, 5), lattice_constants=(1.0, 1.0, 1.0))
>>> sim.set_masses(masses=1.0)
>>> sim.set_random_velocities(temperature=1.0)
>>> sim.set_pair_potential(pair_potential_str='(1-r)**2', r_cut=1.0)
>>> sim.set_neighbor_list(skin=1.0, max_number_of_neighbors=128)
>>> sim.set_integrator(time_step=0.01, target_temperature=1.0, temperature_damping_time=1.0)
>>> print(sim)
Simulation:
    positions: [[0. 0. 0.]] ... [[4. 4. 4.]]
    box_vectors: [5. 5. 5.]
    masses: [[1.]] ... [[1.]]

copy()

Make a deep copy of the simulation

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.positions[:3])
[[0. 0. 0.]
 [0. 0. 1.]
 [0. 0. 2.]]
>>> sim_copy = sim.copy()
>>> sim_copy.positions[0] = [1.0, 0.0, 0.0]
>>> print(sim.positions[:3])
[[0. 0. 0.]
 [0. 0. 1.]
 [0. 0. 2.]]
>>> print(sim_copy.positions[:3])
[[1. 0. 0.]
 [0. 0. 1.]
 [0. 0. 2.]]

from_disk(particle_data='simulation.csv', meta_data='simulation.toml', verbose=False, set_only_particle_data=False) → dict

Load simulation data from disk.

Particle data as CSV file, and meta data as TOML file.

Set simulation box vectors, and particle data. Return meta data from disk as dict.

get_configurational_temperature() → float

Get configurational temperature of the system defined as the ratio of the force squared over the laplacian: T_c=⟨[∇U(r)]²⟩/⟨∇²U(r)⟩ where ⟨…⟩ is average over all particles.

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(f'Configurational temperature: {sim.get_configurational_temperature():.0f}')
Configurational temperature: 0

get_density() → float

Get density of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(f'Density: {sim.get_density():.3f}')
Density: 1.000

get_diameters() → ndarray

Get diameters of particles

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.get_diameters()[:3])
[[1.]
 [1.]
 [1.]]

get_dimensions_of_space() → int

Get dimensions of space

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.get_dimensions_of_space())
3

get_forces() → ndarray

Get forces on all particles

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.positions[0] = [0.5, 0.0, 0.0]  # Shift particle 0 so that the force is not zero
>>> forces = sim.get_forces()
>>> print(f'Force on particle 0: F_x={forces[0, 0]}, F_y={forces[0, 1]}, F_z={forces[0, 2]}')
Force on particle 0: F_x=-1.0, F_y=0.0, F_z=0.0

get_kinetic_energy() → float

Get kinetic energy of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.velocities = np.ones_like(sim.velocities)  # Set all velocities to 1
>>> print(f'Kinetic energy: {sim.get_kinetic_energy():.1f}')
Kinetic energy: 187.5

get_laplacian() → ndarray

Get particle laplacians

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> laplacian = sim.get_laplacian()
>>> print(laplacian[0])  # Laplacian of particle 0
12.0

get_number_of_particles() → int

Return number of particle

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.get_number_of_particles())
125

get_potential_energy() → float

Get total energy of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.scale_box(0.5)
>>> print(sim.get_potential_energy())
167.06442443573454

get_pressure() → float

Get pressure of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_random_velocities(temperature=0.0)
>>> print(sim.get_pressure())
0.0

get_radial_distribution_function(r_bins: ndarray) → tuple[ndarray, ndarray]

Get radial distribution function.

The radial distribution function g(r) is a measure of the probability of finding a particle at a distance r from another particle, relative to the probability of finding a particle in an ideal gas with the same density.

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> r_bins = np.linspace(0.1, 3.0, 100)
>>> g_r = sim.get_radial_distribution_function(r_bins=r_bins)

get_temperature() → float

Get temperature of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(f'Temperature: {sim.get_temperature():.0f}')
Temperature: 1

get_time() → float

Get time of the simulation

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_integrator(time_step=0.01, target_temperature=1.0, temperature_damping_time=0.1)
>>> for _ in range(10):
...     sim.step()
>>> print(sim.get_time())
0.1

get_virial() → float

Get virial of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.get_virial())
0.0

get_volume() → float

Get volume of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.get_volume())
125.0

number_of_particles() → int

Get number of particles

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.number_of_particles())
125

run(steps: int = 1000) → None

Run simulation

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.run(steps=1000)

scale_box(scale_factor: float) → None

Scale box vectors and positions

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.box_vectors)
[5. 5. 5.]
>>> sim.scale_box(0.5)
>>> print(sim.box_vectors)
[2.5 2.5 2.5]

set_density(density: float = 1.0) → None

Set density of the system

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_density(density=0.9)
>>> print(f'Density: {sim.get_density():.3f}')
Density: 0.900

set_integrator(time_step: float = None, target_temperature: float = None, temperature_damping_time: float = None) → None

Set integrator parameters

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_integrator(time_step=0.01, target_temperature=1.0, temperature_damping_time=0.1)

set_masses(masses: float | list | ndarray = 1.0)

Set masses of particles.

The masses can be given as a float, list or ndarray.

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_masses(masses=1.0)
>>> print(sim.masses[:3])
[[1.]
 [1.]
 [1.]]
>>> sim.set_masses(masses=[1.0]*sim.number_of_particles())
>>> print(sim.masses[:3])
[[1.]
 [1.]
 [1.]]
>>> sim.set_masses(masses=np.ones(sim.number_of_particles()))
>>> print(sim.masses[:3])
[[1.]
 [1.]
 [1.]]

set_neighbor_list(skin: float = None, max_number_of_neighbors: int = None, method_str: None | str = None) → None

Update neighbour list

Parameters:

• skin (float) – Skin distance for updating neighbor list.

• max_number_of_neighbors (int) – Maximum number of neighbors for each particle.

• method_str (str) – Method for updating neighbor list.

Examples

>>> from dompap import Simulation
>>> from dompap.neighbor_list import get_number_of_neighbors
>>> sim = Simulation()
>>> sim.set_neighbor_list()
>>> max_number_of_neighbours = np.max(get_number_of_neighbors(sim.neighbor_list))
>>> print(f'Max number of neighbours: {max_number_of_neighbours}')
Max number of neighbours: 26

set_pair_potential(pair_potential_str: str = '(1-r)**2', r_cut: float = 1.0, force_method=None, energy_method=None) → None

Set pair potential.py and force

Parameters:

• pair_potential_str (str) – String representation of the pair potential. Use r as the distance between two particles.

• r_cut (float) – Cutoff distance for the pair potential. I.e the pair potential is zero for r > r_cut.

• force_method (str) – Method for calculating forces. See self._KNOWN_FORCE_METHODS for available methods.

• energy_method (str) – Method for calculating energy. See self._KNOWN_ENERGY_METHODS for available methods.

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_pair_potential(pair_potential_str='(1-r)**2', r_cut=1.0)
>>> sim.set_pair_potential_parameters(sigma=2.0, epsilon=4.0)
>>> print(sim.pair_potential(0.5))
0.25

set_pair_potential_parameters(sigma: float = 1.0, epsilon: float = 1.0) → None

Set potential parameters.

Give sigma and epsilon as floats or callables.

Parameters:

• sigma (float | callable) – Sigma parameter of the pair potential. If callable, it should take two arguments n and m, and return a float.

• epsilon (float | callable) – Epsilon parameter of the pair potential. If callable, it should take two arguments n and m, and return a float.

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_pair_potential_parameters(sigma=1.0, epsilon=1.0)  # Set all sigmas and epsilons to 1.0
>>> print(sim.sigma_func(0, 0))
1.0
>>> print(sim.sigma_func(0, 1))
1.0
>>> def sigma_func(n, m):
...     if n == m:
...         return 1.0
...     else:
...         return 1.2
>>> def epsilon_func(n, m):
...     if n == m:
...         return 1.0
...     else:
...         return 0.5
>>> sim.set_pair_potential_parameters(sigma=sigma_func, epsilon=epsilon_func)
>>> print(sim.sigma_func(0, 0))
1.0
>>> print(sim.sigma_func(0, 1))
1.2
>>> print(sim.epsilon_func(0, 0))
1.0
>>> print(sim.epsilon_func(0, 1))
0.5

set_particle_types(types) → None

Set particle types

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_particle_types(types=1)
>>> print(sim.particle_types[:3])
[[1]
 [1]
 [1]]

set_positions(unit_cell_coordinates: tuple = ([0.0, 0.0, 0.0],), cells: tuple = (5, 5, 5), lattice_constants: tuple = (1.0, 1.0, 1.0))

Set positions of particles

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_positions(unit_cell_coordinates=([0, 0, 0],), cells=(5, 5, 5))
>>> print(sim.positions[:5])
[[0. 0. 0.]
 [0. 0. 1.]
 [0. 0. 2.]
 [0. 0. 3.]
 [0. 0. 4.]]

set_random_velocities(temperature: float = 1.0)

Set velocities from Normal distribution with variance temperature / mass

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.set_random_velocities(temperature=1.0)

step() → None

Make one step in the simulation

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.step()

step_leap_frog() → None

Make one step in the simulation using the NVE Leap-frog integrator

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> sim.step_leap_frog()

to_disk(particle_data='simulation.csv', meta_data='simulation.toml') → None

Save simulation data to disk.

Particle data as CSV file, and meta data as TOML file.

update_neighbor_list(check=True) → None

Update neighbour list if needed

Parameters:

check (bool) – If True, check if neighbour list is old before updating.

wrap_into_box() → None

Wrap all particles into box

Examples

>>> from dompap import Simulation
>>> sim = Simulation()
>>> print(sim.box_vectors)
[5. 5. 5.]
>>> print(sim.image_positions[0])
[0 0 0]
>>> sim.positions[0] = [6.0, 0.0, 0.0]  # Shift particle 0 outside the box
>>> print(sim.positions[0])
[6. 0. 0.]
>>> sim.wrap_into_box()
>>> print(sim.positions[0])
[1. 0. 0.]
>>> print(sim.image_positions[0])
[1 0 0]

dompap.tools

Tools to help with simulations.

dompap.tools.autotune(sim: Simulation, steps=100, test_double_loop=True, smallest_skin=0.1, step_skin=0.1, verbose=False, plot=False) → Simulation

Autotune the neighbor list and force methods for a simulation.

Parameters:

• sim (Simulation) – Simulation to be autotuned.

• steps (int) – Number of steps to run for each skin value.

• test_double_loop (bool) – Whether to test the double loop method.

• smallest_skin (float) – Smallest skin value to test.

• step_skin (float) – Step size for skin values.

• verbose (bool) – Whether to print additional information.

• plot (bool) – Whether to plot results of skin values

Examples

>>> from dompap import Simulation
>>> from dompap.tools import autotune
>>> sim = Simulation()
>>> # Set up a simulation
>>> sim = autotune(sim)