API#

The Simulation Class#

class gamdpy.Simulation(configuration: Configuration, interactions: Interaction | list[Interaction], integrator: Integrator, runtime_actions: list[RuntimeAction], num_timeblocks, steps_per_timeblock, storage: str, compute_plan=None, verbose=False, timing=True, steps_in_kernel_test=1)#

Class for running a simulation.

Parameters:

configuration (gamdpy.Configuration) – The configuration to simulate.
interactions (list of interactions) – Interactions such as pair potentials, bonds, external fields, etc.
integrator (an integrator) – The integrator to use for the simulation.
runtime_actions (list of runtime actions) – Runtime actions such as ScalarSaver, TrajectorySaver or MomentumReset
num_timeblocks (int) – Number of timeblocks to run the simulation. If not 0, then steps_per_timeblock should be set.
steps_per_timeblock (int) – Number of steps per timeblock.
compute_plan (dict) – A dictionary with the compute plan for the simulation. If None, a default compute plan is used.
storage (str) – Storage for the simulation output. Can be ‘memory’ or a filename with extension ‘.h5’.
verbose (bool) – If True, print verbose output.
timing (bool) – If True, timing information is saved.

See also

gamdpy.get_default_sim()

run(verbose=True)#: Run the simulation.

See also

gamdpy.Simulation.timeblocks()

run_timeblocks(num_timeblocks=-1)#

Generator for running the simulation one block at a time.

Parameters:: num_timeblocks (int) – Number of blocks to run. If -1, all blocks are run.

Examples

>>> import gamdpy as gp
>>> sim = gp.get_default_sim()
>>> for block in sim.run_timeblocks(num_timeblocks=3):
...     print(f'{block=}')  # Replace with code to analyze the current configuration
block=0
block=1
block=2

See also

gamdpy.Simulation.run()

status(per_particle=False) → str#

String with the current status Should be executed during the simulation run, see gamdpy.Simulation.timeblocks()

Parameters:: per_particle (bool) – If True, the values are divided by the number of particles in the configuration.
Returns:: A string with the current status of the simulation.
Return type:: str

summary() → str#: Returns a summary string of the simulation run. Should be called after the simulation has been run, see gamdpy.Simulation.timeblocks() or gamdpy.Simulation.run()

The Configuration Class#

class gamdpy.Configuration(D: int, N: int = None, type_names=None, compute_flags=None, ftype=<class 'numpy.float32'>, itype=<class 'numpy.int32'>)#

The configuration class

Store particle vectors (positions, velocities, forces) and scalars (energy, virial, mass …). Also store particle type, image coordinates, and the simulation box.

Parameters:

D (int) – Spatial dimension for the configuration.
N (int [Optional]) – Number of particles. If not set, this will be determined the first time particle data is written to the configuration.

Examples

>>> import gamdpy as gp
>>> conf = gp.Configuration(D=3, N=1000)
>>> print(conf.vector_columns)  # Print names of vector columns
['r', 'v', 'f']
>>> print(conf.scalar_columns) # Print names of scalar columns
['U', 'W', 'K', 'm']
>>> print(conf['r'].shape) # Vectors are stored as (N, D) numpy arrays
(1000, 3)
>>> print(conf['m'].shape) # Scalars are stored as (N,) numpy arrays
(1000,)

Data can be accessed via string keys (similar to dataframes in pandas):

>>> conf['r'] = np.ones((1000, 3))
>>> conf['v'] = 2   # Broadcast by numpy to correct shape
>>> print(conf['r'] + 0.01*conf['v'])
[[1.02 1.02 1.02]
 [1.02 1.02 1.02]
 [1.02 1.02 1.02]
 ...
 [1.02 1.02 1.02]
 [1.02 1.02 1.02]
 [1.02 1.02 1.02]]

A configuration can be specified without setting the number particles, N. In that case N is determined the first time the particle data is written to the configuration:

>>> import numpy as np
>>> conf = gp.Configuration(D=3)
>>> conf['r'] = np.zeros((400, 3))
>>> print(conf['r'].shape)
(400, 3)

copy_to_device()#: Copy all data to device memory

copy_to_host()#: Copy all data to host memory

get_potential_energy() → float#: Get total potential energy of the configuration

get_volume()#: Get volume of simulation box associated with configuration

make_lattice(unit_cell: dict, cells: list, rho: float = None) → None#

Generate a lattice configuration

The lattice is constructed by replicating the unit cell in all directions. Unit cell is a dictonary with fractional_coordinates for particles, and the lattice_constants as a list of unit cell lengths in all directions.

Unit cells are available in gamdpy.unit_cells

Example

>>> import gamdpy as gp
>>> conf = gp.Configuration(D=3)
>>> conf.make_lattice(gp.unit_cells.FCC, cells=[8, 8, 8], rho=1.0)
>>> print(gp.unit_cells.FCC)  # Example of a unit cell dict
{'fractional_coordinates': [[0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]], 'lattice_constants': [1.0, 1.0, 1.0]}

make_positions(N, rho)#

Generate particle positions in D dimensions.

Positions are generated in a simple cubic configuration in D dimensions. Takes the number of particles N and the density rho as inputs

Example

>>> import gamdpy as gp
>>> conf = gp.Configuration(D=3)
>>> conf.make_positions(N=27, rho=0.2)

randomize_velocities(temperature, seed=None, ndofs=None)#: Randomize velocities according to a given temperature. If T <= 0, set all velocities to zero.

save(output: File, group_name: str, mode='w', include_topology=False) → None#

Write a configuration to a HDF5 file

Parameters:

configuration (gamdpy.Configuration) – a gamdpy configuration object
output (h5py.File) – h5 file
group_name (str) – name of the group which will be created in the h5 and in which the configuration will be saved
mode (str) – default value is “w” and corresponds to replacing existing dataset

Example

>>> import os
>>> import h5py
>>> import gamdpy as gp
>>> conf = gp.Configuration(D=3)
>>> conf.make_positions(N=10, rho=1.0)
>>> conf.save(output=h5py.File("final.h5", "w"), group_name="configuration", mode="w")
>>> os.remove("final.h5")       # Removes file (for doctests)
>>> with h5py.File("manyconfs.h5", "a") as fout:
...     conf.save(output=fout, group_name="restarts/restart0000", mode="w")
...     conf.save(output=fout, group_name="restarts/restart0001", mode="w")
...     conf.save(output=fout, group_name="restarts/restart0002", mode="w")
>>> os.remove("manyconfs.h5")       # Removes file (for doctests)

Simulation boxes#

An simulation box object is attached to an configuration object.

class gamdpy.Orthorhombic(D, lengths)#

Standard rectangular simulation box class

Example

>>> import gamdpy as gp
>>> import numpy as np
>>> simbox = gp.Orthorhombic(D=3, lengths=np.array([3,4,5]))

get_dist_sq_dr_function()#: Generates function dist_sq_dr which computes displacement and distance squared for one neighbor

get_dist_sq_function()#: Generates.function dist_sq_function which computes distance squared for one neighbor

get_lengths()#: Return the box lengths as a numpy array

get_volume()#: Return the box volume

get_volume_function()#: Return the box volume

scale(scale_factor)#: Scale the box lengths by scale_factor

class gamdpy.LeesEdwards(D, lengths, box_shift=0.0, box_shift_image=0)#

Simulation box class with LeesEdwards bondary conditions

Example

>>> import gamdpy as gp
>>> import numpy as np
>>> simbox = gp.LeesEdwards(D=3, lengths=np.array([3,4,5]), box_shift=1.0)

get_dist_sq_dr_function()#: Generates function dist_sq_dr which computes displacement and distance for one neighbor

get_dist_sq_function()#: Generates.function dist_sq_function which computes distance squared for one neighbor

get_volume()#: Calculate and return the volume, on the host

get_volume_function()#: Returns the function which calculates the volume of the simulation box

scale(scale_factor)#: Rescale the box by a factor scale_factor, in all directions, if scale_factor is a single float or by a different factor in each direction if scale_factor is an array of length D

Integrators#

class gamdpy.NVE(dt)#

Total energy conserving integrator

Use the Leap-frog algorithm to integrate the equations of motion:

\[ \begin{align}\begin{aligned}v(t+dt/2) &= v(t-dt/2) + f(t) dt / m\\r(t+dt) &= r(t) + v(t+dt/2) dt\end{aligned}\end{align} \]

Parameters:: dt (float) – Time step for the integration.

get_kernel(configuration: Configuration, compute_plan: dict, compute_flags: dict[str, bool], interactions_kernel, verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements performing a number of steps of the integrator

get_params(configuration: Configuration, interactions_params: tuple, verbose=False) → tuple#: Get a tuple with the parameters expected by the associated kernel

class gamdpy.NVT(temperature, tau: float, dt: float)#

The leapfrog algorithm with a Nose-Hoover thermostat

Integrator keeping N (number of particles), V (volume), and T (temperature) constant, using the leapfrog algorithm and a Nose-Hoover thermostat.

Parameters:

temperature (float or function) – The temperature to keep the system at. If a function, it should take time as argument.
tau (float) – The relaxation time of the thermostat.
dt (float) – The time step of the integration.

get_kernel(configuration: Configuration, compute_plan: dict, compute_flags: dict, interactions_kernel, verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements performing a number of steps of the integrator

get_params(configuration: Configuration, interactions_params: tuple, verbose=False) → tuple#: Get a tuple with the parameters expected by the associated kernel

class gamdpy.NVT_Langevin(temperature, alpha: float, dt: float, seed: int)#

NVT Langevin Leap-frog integrator

The langevin thermostat is a stochastic thermostat that keeps the system at a constant temperature:

\[m a = f - \alpha v + \beta\]

where $a$ is the acceleration, $f$ is the force, $v$ is the velocity, $m$ is the mass, $\alpha$ is the friction coefficient, and $\beta$ is a random number drawn from a normal distribution. The implementation uses the leap-frog algorithm described in reference https://arxiv.org/pdf/1303.7011.pdf

Parameters:

temperature (float or function) – Temperature of the thermostat. If a function, it must take a single argument, time, and return a float.
alpha (float) – Friction coefficient of the thermostat.
dt (float) – Time step for the integration.

get_kernel(configuration: Configuration, compute_plan: dict, compute_flags: dict[str, bool], interactions_kernel, verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements performing a number of steps of the integrator

get_params(configuration: Configuration, interactions_params: tuple, verbose=False) → tuple#: Get a tuple with the parameters expected by the associated kernel

class gamdpy.NPT_Atomic(temperature, tau: float, pressure, tau_p: float, dt: float)#

Constant NPT integrator for atomic systems

Integrator keeping N (number of particles), P (pressure), and T (temperature) constant, using the leapfrog algorithm and the thermostat-barostat by G Martyna, DJ Tobias and ML Klein in Molecular Physics 87(5), 1117 (1996), doi:10.1063/1.467468 Note that the thermostat and barostat states defined here are $p_\xi/Q$ and $p_\epsilon/W$ from Eq. 2.9 in the paper.

Parameters:

temperature (float or function) – Target temperature as a function of time or constant
tau (float) – Thermostat relaxation time
pressure (float or function) – Target pressure as a function of time or constant
tau_p (float) – Barostat relaxation time
dt (float) – Timestep size

get_kernel(configuration: Configuration, compute_plan: dict, compute_flags: dict, interactions_kernel, verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements performing a number of steps of the integrator

get_params(configuration: Configuration, interactions_params: tuple, verbose=False) → tuple#: Get a tuple with the parameters expected by the associated kernel

class gamdpy.NPT_Langevin(temperature, pressure, alpha: float, alpha_baro, mass_baro, volume_velocity, barostatModeISO, boxFlucCoord, dt: float, seed: int)#

Constant NPT Langevin integrator NPT Langevin Leap-frog integrator based on N. Grønbech-Jensen and Oded Farago, J. Chem. Phys. 141, 194108 (2014), doi:10.1063/1.4901303

get_kernel(configuration: Configuration, compute_plan: dict, compute_flags: dict[str, bool], interactions_kernel, verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements performing a number of steps of the integrator

get_params(configuration: Configuration, interactions_params: tuple, verbose=False) → tuple#: Get a tuple with the parameters expected by the associated kernel

class gamdpy.SLLOD(shear_rate, dt)#

The SLLOD integrator

Shear an atomic system in the xy-plane using the SLLOD equations.

Parameters:

shear_rate (float) – The shear rate of the system.
dt (float) – The time step of the simulation.

get_kernel(configuration, compute_plan, compute_flags, interactions_kernel, verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements performing a number of steps of the integrator

get_params(configuration, interactions_params, verbose=False)#: Get a tuple with the parameters expected by the associated kernel

class gamdpy.NVU_RT(target_u: float, threshold: float, initial_step: float = 0.1, initial_step_if_high: float = 0.01, step: float = 1, max_steps: int = 20, max_initial_step_corrections: int = 20, max_abs_val: float = 2, eps: float = 1e-07, debug_print: bool = False, mode: Literal['reflection', 'no-inertia', 'reflection-mass_scaling'] = 'reflection-mass_scaling', save_path_u: bool = False, raytracing_method: Literal['parabola', 'parabola-newton', 'bisection'] = 'parabola', float_type: Literal['32', '64'] = '64')#

Potential energy conserving integrator. Calculate the positions by reflecting on the constant potential energy hypersurface and doing Ray Tracing (RT).

Uses parabola approximations, newton-rhapson or bisection to perform raytrcing.

Parameters:

target_u (float) – Target Potential Energy (U_0) to maintain constant along the simulation
threshold (float) – Width of the potential energy “shell” relative to U_0. [Iterative Method] When the potential energy of a certain configuration is within threshold if the potential energy of the fist step, iteration is finished. |U(t) / U_0 - 1| < threshold. It needs to be small enough so that the iterative method is precise enough (it can get very chaotic if it is too high). For example: 1e-6.
initial_step (float, default=0.1) – [Iterative Metod] Initial step in configuration space so that x = positions + d * initial_step is the initial guess for the root algorithm. d is the velocity normalized. For method bisection, it needs to be big enough so that it steps away from the initial position but small enough so that it doesn’t reach out of the surface. If it is two big the algorithm will try to correct it. For method parabola it needs to be as big as possible but it needs to be a point that is below U_0. For example: 0.01. It is better to be based on the density so a good value is something like 0.5 / rho^(1/3).
initial_step_if_high (float, default=0.01) – [Iterative Metod] Initial step if the potential energy of the initial configuration (time == 0) is higher than the target potential energy. It should be high enough to “enter” the potential energy surface U = U_0. For example: same setting as initial step.
step (float, default=1) – [Iterative Method] (only for bisection method) Step to look for a point with u > u0. For example: 1. As in initial it is better for it to based on the density.
max_steps (int, default=20) – [Iterative Method] Maximum calls to the interactions kernel to find a point outside the hypersurface. Good enough value is maybe 10 for parabola method and 100 for bisection method.
max_initial_step_corrections (int, default=20) – [Iterative Method] If initial_step is too big the algorithm will try to correct it this amount of times. At the nth correction step initial_step_n = initial_step_0 * (1/2)^n. In the first iteration if U > U0 then initial step is corrected to make it bigger (s_n = s_0 * 2^n). For example: 10 (max correction will be approximately 1000).
max_abs_val (float, default=2) – [Iterative Metod] (only for bisection) Some potential energy functions increas rapidly at a certain configurations. To prevent numerical errors, if the absolute potential energy of a given configuration, |U|, is above U_0 * max_abs_val, the iterative method will discard that configuration as invalid and reach “less far” into configurational space. For example: 2.
eps (float, default=1e-7) – [Iterative Method] Because of numerical inaccuracies, it could happen that the same value calculated twice once is positive and once is negative. Values of the potential energy relative to the target potential enery with |x| < eps are considered neither positive or negative in the algo. Needs to be smaller than threshold.
debug_print (bool, default=False) – If a root is not found, the 0th thread prints useful debugging information if debug_print is enabled.
mode ({"reflection", "no-inertia", "reflection-mass_scaling"}, default = "reflection-mass_scaling") – Mode to perform the reflection. reflection-mass_scaling applies a correction that takes into account the mass of each particle. If the setup does not include particles with different masses, reflection will be faster (not that much). no-inertia is a testing feature: instead of reflecitng velocities in the hyper surface the new velocities follow the direction of the normal vector (the force).
save_path_u (optional, default=False) – Save the potential energy between two consecutive points in the iteration
raytracing_method ({"parabola", "parabola-newton", "bisection"}, default = "parabola") – Methodology to find the next point in the surface with same potential energy
float_type ({"64", "32"}, default = "64") – Float type for the potential energy. Higher threshold can work with float32

get_kernel(configuration, compute_plan, compute_flags, interactions_kernel, verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements performing a number of steps of the integrator

get_params(configuration, interaction_params, verbose=False)#: Get a tuple with the parameters expected by the associated kernel

Interactions#

Pair potentials#

class gamdpy.PairPotential(pairpotential_function, params, max_num_nbs, exclusions=None)#

Pair potential

check_datastructure_validity() → bool#: Interactions which have an internal data structure (think: PairPotential and its NbList) should overwrite this method with one that checks the validity of it, and throws an error if not valid (Later: repair it and signal rerun of timeblock)

get_kernel(configuration: Configuration, compute_plan: dict, compute_flags: dict[str, bool], verbose=False)#: Get a kernel (or python function depending on compute_plan[“gridsync”]) that implements calculation of the interaction

get_params(configuration: Configuration, compute_plan: dict, verbose=False) → tuple#: Get a tuple with the parameters expected by the associated kernel

Functions (pair potentials)#

gamdpy.LJ_12_6(dist, params)#

The 12-6 Lennard-Jones potential

See gamdpy.apply_shifted_potential_cutoff() for a usage example.

\[ \begin{align}\begin{aligned}u(r) = A_{12} r^{-12} + A_6 r^{-6}\\u'(r) = -12*A_{12} r^{-13} - 6*A_6 r^{-7}\end{aligned}\end{align} \]

Parameters:

dist (float) – Distance between particles
params (array-like) – A₁₂, A₆

Returns:

u (float) – Potential energy
s (float) – Force multiplier, -u’(r)/r
umm (float) – Second derivative of potential energy

gamdpy.LJ_12_6_sigma_epsilon(dist, params)#

The 12-6 Lennard-Jones potential

\[u(r) = 4\epsilon( (r/\sigma)^{-12} - (r/\sigma)^{-6} )\]

This is the same as the gamdpy.LJ_12_6() potential, but with $\sigma$ (sigma) and $\epsilon$ (epsilon) as parameters.

Parameters:

dist (float) – Distance between particles
params (array-like) – σ, ε

gamdpy.harmonic_repulsion(dist, params)#

The harmonic repulsion pair potential

\[u(r) = ½\epsilon(1-r/\sigma)^2\]

for $r<\sigma$ and zero otherwise. Parameters: ε=epsilon, σ=cut. Note that this potential is naturally truncated at r=σ.

Parameters:

dist (float) – Distance between particles
params (array-like) – ε, σ

gamdpy.hertzian(dist, params)#

Hertzian potential

\[u(r) = \epsilon(1-r/\sigma)^\alpha\]

for $r<\sigma$ and zero otherwise. Parameters: ε=epsilon, α=alpha, σ=cut. Note that this potential is naturally truncated at r=σ. For Hertzian disks chose α=7/2, and for Hertzian spheres chose α=5/2

Parameters:

dist (float) – Distance between particles
params (array-like) – ε, α, σ

gamdpy.SAAP(dist, params)#

The SAAP potential is a pair potential for noble elements, its parameters are derived from fitting on data obtained from ab initio methods (coupled cluster):

\[u(x=r/\sigma)/\epsilon = (a_0\exp(a_1 x)/x + a_2\exp(a_3 x) + a_4) / (1+a_5 x^6)\]

The six $a_i$ parameters are given in units of eps and sigma.

Reference: https://doi.org/10.1063/1.5085420

Parameters:

dist (float) – Distance between particles
params (array-like) – a₀, a₁, a₂, a₃, a₄, a₅, σ, ε

Functions (bonds)#

gamdpy.harmonic_bond_function(dist: float, params: ndarray) → tuple#

Harmonic bond potential

\[u(r) = \frac{1}{2} k (r - r_0)^2\]

Parameters:

dist (float) – Distance between particles
params (array-like) – r₀, k

Returns:

u (float) – Potential energy
s (float) – Force multiplier, -u’(r)/r
umm (float) – Second derivative of potential energy

See also

gamdpy.Bonds

Generators#

Generators return a function that can be used to calculate the potential energy and the force between two particles.

gamdpy.add_potential_functions(potential_1, potential_2)#

Add two potential functions into a single potential function

Note that the two potential functions will have the same cut-off, by convention always stored as the last entry in params. The potential_1 cannot explicitly depend on cut-off (last entry in params). The potential_2 should know where to look for its first parameter in the list of parmeters (params), see e.g. first_parameter in gamdpy.make_IPL_n().

Parameters:

potential_1 (callable) – a function that calculates a pair-potential: u, s, umm = potential_1(dist, params)
potential_2 (callable) – a function that calculates a pair-potential: u, s, umm = potential_2(dist, params)

Returns:

potential – a function implementing the sum of potential_1 and potential_2

Return type:

callable

Example

Below we make the 12-6 Lennard-Jones potential by adding two inverse power-law potentials.

>>> import gamdpy as gp
>>> LJ = gp.add_potential_functions(gp.make_IPL_n(12), gp.make_IPL_n(6, first_parameter=1))
>>> params = A12, A6, cut = 4.0, -4.0, 2.5
>>> dist = 2**(1/6) # Minima of LJ potential
>>> LJ(dist,params)[0] # Pair energy in minima
-1.0

gamdpy.make_potential_function_from_sympy(ufunc, param_names) → callable#

Make a potential function from a sympy expression

Known to result in slow code. Use for testing and prototyping.

Parameters:

ufunc (sympy expression) – The potential energy expression in Sympy’s symbolic form It has to be a radial function and the pair distance symbol shuould be r
param_names (list) – List of parameters

Returns:

potential_function – A function that calculates the potential energy, force multiplier, and second derivative of the potential energy u, s, umm = potential_function(dist, params)

Return type:

callable

gamdpy.make_LJ_m_n(m: float, n: float) → callable#

Mie Potential

Also known as the generalized Lennard-Jones potential:

\[u(r) = A_m r^{-m} - A_n r^{-n}\]

Returns:: potential_function – A function that calculates the Mie potential, u, s, umm = potential_function(dist, params). where params = [A_m, A_n]
Return type:: callable

gamdpy.make_IPL_n(n: float, first_parameter: int = 0) → callable#

Inverse Power Law Potential

\[u(r) = A_n r^{-n}\]

Parameters:

n (float) – Exponent in the potential
first_parameter (int) – The index of the first parameter in the list of parameters. See usage in gamdpy.add_potential_functions().

Returns:

potential_function – A function that calculates the IPL potential, u, s, umm = potential_function(dist, params). where params = [A_n]

Return type:

callable

Modifies#

Modifies are typically used to smoothly truncate the potential at a certain distance.

gamdpy.apply_shifted_potential_cutoff(pair_potential: callable) → callable#

Apply shifted potential cutoff to a pair-potential function

If the input pair potential is $u(r)$, then the shifted potential is $u(r) - u(r_{c})$, where $r_c$ is the cutoff distance.

Note: calls original potential function twice, avoiding changes to params.

Parameters:: pair_potential (callable) – A function that calculates a pair-potential: u, s, umm = pair_potential(dist, params)
Returns:: pair_potential – A function where shifted_potential_cutoff is applied to original function
Return type:: callable

Example

The following example demonstrates how to use this function to set up the Lenard-Jones 12-6 potential truncated and shifted to zero at the cutoff distance of 2.5:

>>> import gamdpy as gp
>>> pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6)
>>> A12, A6, cut = 1.0, 1.0, 2.5
>>> pair_pot = gp.PairPotential(pair_func, params=[A12, A6, cut], max_num_nbs=1000)

gamdpy.apply_shifted_force_cutoff(pair_potential)#

Apply shifted force cutoff to a pair-potential function

If the input pair potential is $u(r)$, then the shifted force potential is $u(r) - u(r_{c}) + s(r_{c})(r - r_{c})$, where $r_c$ is the cutoff distance, and $s(r) = -u'(r)/r$.

Note: calls original potential function twice, avoiding changes to params

Parameters:: pair_potential (callable) – a function that calculates a pair-potential: u, s, umm = pair_potential(dist, params)
Returns:: potential – a function where shifted force cutoff is applied to original function
Return type:: callable

Fixed interactions#

Classes#

class gamdpy.Bonds(bond_potential, potential_params, indices)#

Fixed bond interactions between particles, such as harmonic bonds or FENE bonds.

Parameters:

bond_potential (function) – A function that takes the distance between two particles and the bond type as arguments and returns the potential energy, force and laplacian. See :func:gamdpy.potential_functions.harmonic_bond_function for an example.
potential_params (list) – A list of parameters for each bond type. Each entry is a list of parameters for a specific bond type.
indices (list) – A list of lists, each containing the indices of the two particles involved in a bond and the bond

See also

gamdpy.harmonic_bond_function: Harmonic bond potential

class gamdpy.Angles(potential, indices, parameters)#

class gamdpy.Tether#

Connect particles to anchor-points in space with a harmonic spring force.

Parameters:

either (Points and spring constants are defined using)
constants ((b) A list of particle types to be tethered and associated list of spring)
constants
examples/tethered_particles.py (See)

class gamdpy.Gravity#

Adding a gravitational-like force on particles.

Parameters:

forces ((b) A list of particle types one which the forces act and associated list of)
forces
x-direction (At the moment the force will act in the)
examples/poiseuille.py (See)

class gamdpy.Relaxtemp#

Thermostat using simple kinetic temperature relaxation of each particles.

Parameters:

time ((a) A list of particle indices to be thermostated and associated list of relaxation)
times ((b) A list of particle types to be thermostated and associated list of relaxation)
examples/poiseuille.py (See)

Generators#

gamdpy.make_planar_calculator(configuration, potential_function) → callable#

Returns a function that calculates a planar interaction for particles

This function is used to create a planar interaction such as a smooth wall, gravity or an electric field.

gamdpy.setup_planar_interactions(configuration, potential_function, potential_params_list, particles_list, point_list, normal_vector_list, compute_plan, verbose=True) → dict#: Returns a dictionary with planar interactions

gamdpy.make_fixed_interactions(configuration, fixed_potential, compute_plan, verbose=True)#: Generate a kernel for fixed interactions between particles.

Runtime Actions#

class gamdpy.TrajectorySaver(include_simbox=False, verbose=False, compression='gzip', compression_opts=4)#: Runtime action for saving configurations during timeblock. Does logarithmic saving.

class gamdpy.ScalarSaver(steps_between_output: int = 16, compute_flags=None, verbose=False, compression='gzip', compression_opts=4)#: Runtime action for saving scalar data (such as thermodynamic properties) during a timeblock every steps_between_output time steps.

class gamdpy.RestartSaver(timeblocks_between_restarts=1)#: Runtime action for saving restarts, ie. the current configuration at beginning of every timebock .

class gamdpy.MomentumReset(steps_between_reset: int)#: Runtime action that sets the total momentum of configuration to zero every steps_between_reset time step.

class gamdpy.StressSaver(steps_between_output: int = 16, compute_flags=None, verbose=False)#: Runtime action for saving stress tensor(a D x D matrix) during a timeblock every steps_between_output time steps.

Calculators#

class gamdpy.CalculatorRadialDistribution(configuration, bins, compute_plan=None, ptype=None)#

Calculator class for the radial distribution function, g(r)

Parameters:

configuration (gamdpy.Configuration) – The configuration object for which the radial distribution function is calculated.
bins (int) – The number of bins in the radial distribution function.
compute_plan (dict)
ptype (optional) – If specified, array ptypes used to calculate g(r). If not specfified default is configuration.ptype

Example

>>> import gamdpy as gp
>>> sim = gp.get_default_sim()
>>> calc_rdf = gp.CalculatorRadialDistribution(sim.configuration, bins=1000)
>>> for _ in sim.run_timeblocks():
...     calc_rdf.update()      # Current configuration to rdf
>>> rdf_data = calc_rdf.read() # Read the rdf data as a dictionary
>>> r = rdf_data['distances']  # Pair distances
>>> rdf = rdf_data['rdf']      # Radial distribution function

read()#

Read the radial distribution function

Returns:: A dictionary containing the distances and the radial distribution function.
Return type:: dict

save_average(output_filename='rdf.dat', save_ptype=False) → None#

Save the average radial distribution function to a file

Parameters:

output_filename (str) – The name of the file to which the radial distribution function is saved.
save_ptype (bool) – Save the type of each particle in a file name ptype_* (default False)

update()#: Update the radial distribution function with the current configuration.

class gamdpy.CalculatorStructureFactor(configuration: Configuration, n_vectors: ndarray = None, atomic_form_factors: ndarray = None, backend='CPU multi core')#

Calculator class for the static structure factor, S(q). The calculation is done for several ${\bf q}$ vectors given by

\[{\bf q} = (2\pi n_x/L_x, 2\pi n_y/L_y, ...)\]

where $n=(n_x, n_y, ...)$ is a D-dimensional vector of integers and $L_x$, $L_y$, … are the box lengths in the $x$, $y$, … directions, respectively. Note that box length are assumed to be constant during the simulation (as in a NVT simulation). The collective density $\rho_{\bf q}$ is calculated as

\[\rho_{\bf q} = \frac{1}{\sqrt{N}} \sum_{n} f_n \exp(-i {\bf q}\cdot {\bf r}_n)\]

where $x_n$ is the position of particle $n$, and $f_n$ is the atomic form factor for that particle (one by default). The normalization constant is

\[N = \sum_{n} f_n.\]

From this, the static structure factor is defined as

\[S({\bf q}) = |\rho_{\bf q}|^2.\]

The method update() updates the structure factor with the current configuration. The method read() returns the structure factor for the q vectors in the q_direction.

Parameters:

configuration (gamdpy.Configuration) – The configuration object to calculate the structure factor for.
q_max (float or None) – The maximum value of the q vectors.
n_vectors (numpy.ndarray or None) – n-vectors defining q-vectors. The shape of n_vectors, if specified, must be (N, D) where N is the number of q vectors and D is the number of dimensions. If None, then use generate_q_vectors method.
atomic_form_factors (numpy.ndarray or None) – The atomic form factors, $f_n$. If None (default), then the atomic form factors are set to 1. Can be given as an array of floats, one for each atom.
backend (str) – The backend to use for the calculation. Either ‘CPU multi core’ or ‘CPU single core’.

See also

CalculatorRadialDistribution

generate_q_vectors(q_max: float)#: Generate q-vectors inside a sphere of radius q_max

read(bins) → dict#

Return the structure factor S(q) for the q vectors in the q_direction.

Parameters:: bins (int | None) – If bins is an integer, the data is binned (ready to be plotted). If bins is None, the raw S(q) data is returned.
Returns:: A dictionary containing the q vectors, the q lengths, the structure factor S(q), the collective density rho_q, and the number of q vectors in each bin. Output depends on the value of the bins parameter.
Return type:: dict

save_average(bins=None, output_filename='sq.dat')#: Save average structure factors to a file.

update() → None#: Update the structure factor with the current configuration.

class gamdpy.CalculatorWidomInsertion(configuration, pair_potential, temperature, ghost_positions, ptype_ghost=0, compute_plan=None, backend='GPU')#

Calculator class for Widom’s particle insertion method for computing chemical potentials

This class is used to calculate the chemical potential of a not to dense fluid using Widom’s particle insertion method. The excess chemical potential, μᵉˣ, fluid can be computed as μᵉˣ = -kT ln 〈exp(-ΔU/kT)〉 where ΔU is the energy difference between the system with and without a ghost particle. Here, 〈…〉 is an average for all possible positions of the ghost particle (sometimes writte as an integral over space).

The method should be used on a NVT trajectory in equlibrium.

Parameters:

configuration (gamdpy.Configuration) – The configuration object representing the system.
pair_potential (gamdpy.PairPotential) – The pair potential object representing the interaction between particles.
temperature (float) – The temperature of the system.
ghost_positions (ndarray) – The positions of the ghost particles for which the chemical potential is to be calculated.
ptype_ghost (int, optional) – The particle type of the ghost particles. Default is 0.
compute_plan (dict, optional) – The compute plan for the system. Default is None (then a default plan is used).
backend (str, optional) – The backend to use for the calculations. Default is ‘GPU’. Supported backends are ‘CPU’ (testing) and ‘GPU’ (recomended).

Example

>>> import numpy as np
>>> import gamdpy as gp
>>> sim = gp.get_default_sim()  # Replace with your own equbriliated simulation
>>> pair_pot = sim.interactions[0]
>>> num_ghost_particles = 500_000
>>> ghost_positions = np.random.rand(num_ghost_particles, sim.configuration.D) * sim.configuration.simbox.get_lengths()
>>> calc_widom = gp.CalculatorWidomInsertion(sim.configuration, pair_pot, sim.integrator.temperature, ghost_positions)
>>> for block in sim.run_timeblocks():
...     calc_widom.update()
>>> calc_widom_data = calc_widom.read()  # Dictionary with chemical potential and more
>>> calc_widom_data.keys()
dict_keys(['chemical_potential', 'boltzmann_factors', 'chemical_potentials'])
>>> mu_ex = calc_widom_data['chemical_potential']  # Estimated excess chemical potential

read()#

Read data

Return the current chemical potential, average Boltzmann factors, and timeblock-specific chemical potentials for the ghost particles.

Returns:

A dictionary containing the following keys:

’chemical_potential’: float
The best estimate of the overall chemical potential for the system.
’boltzmann_factors’: ndarray
The average Boltzmann factors for each ghost particle, representing the statistical weights based on their interaction energies.
’chemical_potentials’: ndarray
An array of estimated chemical potentials for each timeblock, providing insight into the variability of the chemical potential over time.

Return type:

dict

update()#: Update state the current configuration.

class gamdpy.CalculatorHydrodynamicCorrelations(configuration, dtsample, ptype=0, nwaves=10, lvec=100, verbose=False)#

Calculates the tranverse current autocorrelation function (jacf) and the longitudinal density autocorrelation function (dacf).

Example: See hydrocorr.py

Initialization variables - configuration: Instance of the configuration class - dtsample: Time between sampling (= integrator time step * steps_per_timeblock) - ptype: The particle type for which the correlations are calculated (default 0) - nwaves: Number of wavevectors (default 10) - lvec: Length of the time vector - sample time span then dtsample*lvec (defulat 100) - verbose: Write a bit to screen (default False)

Output:: With the method read() data is written to files and returned to user

class gamdpy.CalculatorHydrodynamicProfile(configuration, ptype, bins=100, profdir=2, veldir=0, verbose=True)#

Calculates the density, streaming velocity, and kinetic temperature profiles.

Example: See atomistic_wall example

Initialization variables: - configuration: Instance of the configuration class - ptype: The particle type for which the profile is calculated - bins: Number of bins to use in the profiles; more bins higher resolution. (default 100) - profdir: Profile spatial direction (default 2 - or z - direction) - veldir: The streaming velocity component used (default 0 - or x - direction) - verbose: If true print some information (default True)

Output:: With the method read() you can retrieve the current data. By default a data file is printed with the same

Tools and helper functions#

Input and Output#

The TrajectoryIO class#

class gamdpy.tools.TrajectoryIO(name='')#

This class handles the loading and saving of simulation data. When the class can be instantiated with an output from a previous simulation. The data can be saved in self.h5 in the same format of the sim.output object. When the class is instantiated without an input, self.h5 is None and can be assigned afterward (used to save output from memory simulation)

Parameters:: name (str) – Name of the file or folder to read output from. Can be a gamdpy .h5 output or a rumd3 TrajectoryFiles folder.

Examples

>>> import gamdpy as gp
>>> import h5py
>>> output = gp.tools.TrajectoryIO("examples/Data/LJ_r0.973_T0.70_toread.h5").get_h5()
Found .h5 file (examples/Data/LJ_r0.973_T0.70_toread.h5), loading to gamdpy as output dictionary
>>> nblocks, nconfs, N, D = output['trajectory_saver/positions'].shape
>>> print(f"Output file examples/Data/LJ_r0.973_T0.70.h5 contains a simulation of {N} particles in {D} dimensions")
Output file examples/Data/LJ_r0.973_T0.70.h5 contains a simulation of 2048 particles in 3 dimensions
>>> print(f"The simulation output is divided into {nblocks} blocks, each of them with {nconfs} configurations")
The simulation output is divided into 32 blocks, each of them with 12 configurations
>>> output = gp.tools.TrajectoryIO("examples/Data/NVT_N4000_T2.0_rho1.2_KABLJ_rumd3/TrajectoryFiles").get_h5()   # Read from rumd3
Found rumd3 TrajectoryFiles, loading to rumpdy as output dictionary
>>> nblocks, nconfs, N, D = output['trajectory_saver/positions'].shape
>>> print(f"File examples/Data/NVT_N4000_T2.0_rho1.2_KABLJ_rumd3/TrajectoryFiles contains a simulation of {N} particles in {D} dimensions")
File examples/Data/NVT_N4000_T2.0_rho1.2_KABLJ_rumd3/TrajectoryFiles contains a simulation of 4000 particles in 3 dimensions
>>> print(f"The simulation output is divided into {nblocks} blocks, each of them with {nconfs} configurations")
The simulation output is divided into 2 blocks, each of them with 24 configurations

get_h5() → File#: Returns self.h5

load_h5(name: str) → File#: Makes self.h5 a view of the .h5 files

load_rumd3(name: str) → File#: Reads a rumd3 TrajectoryFiles folder and convert it into gamdpy .h5 output. Assumes RectangularSimulationBox, corresponding to gamdpy’s Orthorhombic. This function returns a memory .h5

save_h5(name: str)#: This method saves self.h5 to disk. It can be used to save sim.output to file if class is initialized without arguments and then self.h5 = sim.output . By default output is compressed. This can be avoided setting self.compression_type = “gzip” and self.compression_opts=0 .

IO functions#

gamdpy.tools.save_configuration(configuration, filename: str, format='xyz', append=False)#

Saves configuration to file.

Helpful for, e.g., VMD visualization This current version only supports the standard xyz-format.

Example

>>> import os
>>> import gamdpy as gp
>>> import numpy as np
>>> conf = gp.Configuration(D=3, N=10)
>>> gp.tools.save_configuration(configuration=conf, filename="final.xyz", format="xyz")
>>> os.remove("final.xyz")      # Removes file (for doctests)

gamdpy.configuration_to_hdf5(configuration: Configuration, filename: str, meta_data=None) → None#

Write a configuration to a HDF5 file

Parameters:

configuration (gamdpy.Configuration) – a gamdpy configuration object
filename (str) – filename of the output file .h5
meta_data (str) – not used in the function so far (default None)

Example

>>> import os
>>> import gamdpy as gp
>>> conf = gp.Configuration(D=3)
>>> conf.make_positions(N=10, rho=1.0)
>>> gp.configuration_to_hdf5(configuration=conf, filename="final.h5")
>>> os.remove("final.h5")       # Removes file (for doctests)

gamdpy.configuration_from_hdf5(filename: str, reset_images=False, compute_flags=None) → Configuration#

Read a configuration from a HDF5 file

Parameters:

filename (str) – filename of the input file .h5
reset_images (bool) – if True set the images to zero (default False)

Returns:

configuration – a gamdpy configuration object

Return type:

gamdpy.Configuration

Example

>>> import gamdpy as gp
>>> conf = gp.configuration_from_hdf5("examples/Data/final.h5")
>>> print(conf.D, conf.N, conf['r'][0])     # Print number of dimensions D, number of particles N and position of first particle
3 10 [-0.7181449 -1.3644753 -1.5799187]

gamdpy.configuration_to_rumd3(configuration: Configuration, filename: str) → None#

Write a configuration to a RUMD3 file

Parameters:

configuration (gamdpy.Configuration) – a gamdpy configuration object
filename (str) – filename of the output file .xyz.gz

Example

>>> import os
>>> import gamdpy as gp
>>> conf = gp.Configuration(D=3)
>>> conf.make_positions(N=10, rho=1.0)
>>> gp.configuration_to_rumd3(configuration=conf, filename="restart.xyz.gz")
>>> os.remove("restart.xyz.gz")       # Removes file (for doctests)

gamdpy.configuration_from_rumd3(filename: str, reset_images=False, compute_flags=None) → Configuration#

Read a configuration from a RUMD3 file

Parameters:: filename (str) – filename of the output file .xyz.gz
Returns:: configuration – a gamdpy configuration object
Return type:: gamdpy.Configuration

Example

>>> import gamdpy as gp
>>> conf = gp.configuration_from_rumd3("examples/Data/NVT_N4000_T2.0_rho1.2_KABLJ_rumd3/TrajectoryFiles/restart0000.xyz.gz")
>>> print(conf.D, conf.N, conf['r'][0])     # Print number of dimensions D, number of particles N and position of first particle
3 4000 [ 7.197245   6.610052  -4.7467813]

gamdpy.configuration_to_lammps(configuration, timestep=0) → str#

Convert a configuration to a string formatted as LAMMPS dump file

Parameters:

configuration (gamdpy.Configuration) – a gamdpy configuration object
timestep (float) – time at which the configuration is saved

Returns:

string formatted as LAMMPS dump file

Return type:

str

Example

>>> import gamdpy as gp
>>> conf = gp.Configuration(D=3)
>>> conf.make_positions(N=10, rho=1.0)
>>> lmp_dump = gp.configuration_to_lammps(configuration=conf)

Post-analysis tools#

gamdpy.tools.calc_dynamics(trajectory, first_block, qvalues=None)#

Compute dynamical properties from a saved trajectory HDF5 file.

This function processes blocks of saved configurations to evaluate time‐dependent dynamical observables, including the mean square displacement (MSD), the non‐Gaussian parameter (alpha2), and the self‐intermediate scattering function (Fs), for one or more particle types.

Parameters:

trajectory (h5py.File object in the gamdpy style)
first_block (int) – Index of the first block to use as the reference origin.
qvalues (float or array‐like of shape (num_types,), optional) – Wavevector magnitudes at which to compute the self‐intermediate scattering function Fs. If a single float is provided, it is broadcast to all particle types. If None, Fs is not computed (remains zero).

Returns:

results – Dictionary containing dynamcal data.

Return type:

dict

Examples

For command‐line usage, see:: $ python -m gamdpy.tools.calc_dynamics -h

Usage within a Python script:

>>> import gamdpy as gp
>>> sim = gp.get_default_sim()  # Replace with your simulation object
>>> for block in sim.run_timeblocks(): pass
>>> dynamics = gp.calc_dynamics(sim.output, first_block=0, qvalues=7.5)
>>> dynamics.keys()
dict_keys(['times', 'msd', 'alpha2', 'qvalues', 'Fs', 'count'])

Mathematical functions#

The below returns functions that can be executed fast in a GPU kernel. As an example, they can be used to set a time-dependent target temperature.

gamdpy.make_function_constant(value)#: Return a function that returns a constant value

gamdpy.make_function_ramp(value0, x0, value1, x1)#

Return a function that ramps linearly between two values

Return a function that returns a constant value for x<x0, linearly ramps from value0 to value1 for x0<=x<=x1, and returns value1 for x>x1.

gamdpy.make_function_sin(period, amplitude, offset)#: Return a function that returns a sin function with given period, amplitude and offset

Extract data#

gamdpy.extract_scalars(data, column_list, first_block=0, D=3)#

Extracts scalar data from simulation output.

Parameters:

data (dict) – Output from a Simulation object.
column_list (list of str)
first_block (int) – Index of the first timeblock to extract data from.
D (int) – Dimension of the simulation.

Returns:

Tuple of 1D numpy arrays containing the extracted scalar data.

Return type:

tuple

Example

>>> import numpy as np
>>> import gamdpy as gp
>>> sim = gp.get_default_sim()  # Replace with your simulation object
>>> for block in sim.run_timeblocks(): pass
>>> U, W = gp.extract_scalars(sim.output, ['U', 'W'], first_block=1)

Miscellaneous#

gamdpy.select_gpu()#: Select the first available GPU (necessary on some multi-gpu systems).

gamdpy.get_default_sim()#

Return a sim object of the single component LJ crystal in the NVT ensemble. The purpose of this function is to provide a default simulation for testing and simplifying examples.

Example

>>> import os
>>> os.environ['NUMBA_CUDA_LOW_OCCUPANCY_WARNINGS'] = '0'   # Removes warnings from low occupacy (optional)
>>> import gamdpy as gp
>>> sim = gp.get_default_sim()

gamdpy.get_default_compute_plan(configuration)#

Return a default compute_plan The default compute_plan is a dictionary with a set of parameters specifying how computations are done on the GPU. The returned plan depends on the number of particles, and properties of the GPU. The keys of the dictionary are:

‘pb’: particle per thread block
‘tp’: threads per particle
‘gridsync’: Boolean indicating if syncronization should be done by grid.sync() calls
‘skin’: used when updating nblist
‘UtilizeNIII’: Boolean indicating if Newton’s third law (NIII) should be utilized (see pairpotential_calculator).
‘nblist’ : ‘N squared’ (default) or ‘linked lists’. Determines algorithm updating nblist

gamdpy.get_default_compute_flags()#

gamdpy.plot_molecule(top, positions, particle_types, filename='molecule.pdf', block=False)#

This function write a pdf file with a drawing of the molecule

Parameters:

top (gamdpy topology object)
positions (list or numpy array with positions of all atoms)
particle_types (types of the molecule)
filename (name of the output pdf file, default is molecule.pdf)
block (boolean, default False. If True shows plot and blocks script until display window is closed)

gamdpy.tools.print_h5_structure(node, indent=0)#

Recursively print groups and datasets with metadata of an h5 file.

Example

>>> import gamdpy as gp
>>> sim = gp.get_default_sim()
>>> for _ in sim.run_timeblocks(): pass
>>> gp.tools.print_h5_structure(sim.output)
initial_configuration/ (Group)
    ptype  (Dataset, shape=(2048,), dtype=int32)
    r_im  (Dataset, shape=(2048, 3), dtype=int32)
    scalars  (Dataset, shape=(2048, 4), dtype=float32)
    topology/ (Group)
        angles  (Dataset, shape=(0,), dtype=int32)
        bonds  (Dataset, shape=(0,), dtype=int32)
        dihedrals  (Dataset, shape=(0,), dtype=int32)
        molecules/ (Group)
    vectors  (Dataset, shape=(3, 2048, 3), dtype=float32)
scalar_saver/ (Group)
    scalars  (Dataset, shape=(8, 64, 3), dtype=float32)
trajectory_saver/ (Group)
    images  (Dataset, shape=(8, 12, 2048, 3), dtype=int32)
    positions  (Dataset, shape=(8, 12, 2048, 3), dtype=float32)

gamdpy.tools.print_h5_attributes(obj, path='/')#

Recursively print attrs of every group/dataset of an h5 file.

Example

>>> import gamdpy as gp
>>> sim = gp.get_default_sim()
>>> for _ in sim.run_timeblocks(): pass
>>> gp.tools.print_h5_attributes(sim.output)
Attributes at /:
    - dt: 0.005
    - script_content: ...
    - script_name: ...
Attributes at /initial_configuration/:
    - simbox_data: [12.815602 12.815602 12.815602]
    - simbox_name: Orthorhombic
Attributes at /initial_configuration/scalars:
    - scalar_columns: ['U' 'W' 'K' 'm']
Attributes at /initial_configuration/topology/molecules/:
    - names: []
Attributes at /initial_configuration/vectors:
    - vector_columns: ['r' 'v' 'f']
Attributes at /scalar_saver/:
    - compression_info: gzip with opts 4
    - scalar_names: ['U' 'W' 'K']
    - steps_between_output: 16
Attributes at /trajectory_saver/:
    - compression_info: gzip with opts 4

API#

The Simulation Class#

The Configuration Class#

Simulation boxes#

Integrators#

Interactions#

Pair potentials#

Functions (pair potentials)#

Functions (bonds)#

Generators#

Modifies#

Fixed interactions#

Classes#

Generators#

Runtime Actions#

Calculators#

Tools and helper functions#

Input and Output#

The TrajectoryIO class#

IO functions#

Post-analysis tools#

Mathematical functions#

Extract data#

Miscellaneous#

This Page