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, timing=True, steps_in_kernel_test=1)

Class for running a simulation.

Parameters:

• configuration (Configuration) – The configuration to simulate.

• interactions (one or a list of interactions) – One or a list of interactions such as pair potentials, bonds, external fields, etc. See the Interactions section for more.

• integrator (an integrator) – The integrator to use for the simulation. See the Integrators section for more.

• runtime_actions (list of runtime actions) – List of runtime actions. See the Runtime Actions section for more.

• num_timeblocks (int) – Number of timeblocks of run the simulation.

• steps_per_timeblock (int) – Number of steps in each 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’.

• timing (bool) – If True, timing information is saved.

See also

gamdpy.get_default_sim()

autotune(time_steps=None, parameters_to_optimize=None, include_linked_lists=True)

Autotune the simulation parameters for most efficient calculations on the current machine

compress(desired_rho, steps_per_rescale, relative_change)

Compress simulation to the density ‘desired_rho’. This is done by a sequence of small adjustments of the simulation box, each followed by a short simulation.

Parameters:

• desired_rho (float) – The final density to scale the simulation to.

• steps_per_rescale (int) – The number of timesteps performed after each adjustments of the simulation box Can not be larger than the ‘steps_per_timeblock’ that the simulation object was created with

• relative_change (float) – The relative increase in density for each adjustments of the simulation box. The last adjustment is computed to hit the desired density exactly

run(verbose=True) → None

Run all blocks of the simulation. See gamdpy.Simulation.run_timeblocks() for an open loop version.

run_timeblocks(num_timeblocks=-1)

Generator for running the simulation one block at a time. The state of the simulation object is updated between block, and data is copied to the host for (optional) data analysis.

Parameters:

num_timeblocks (int) – Number of blocks to run. All blocks are run if -1 (recommended).

Yields:

int – Index of the current block.

Examples

>>> import gamdpy as gp
>>> sim = gp.get_default_sim()  # Replace this with your own simulation object
>>> for block in sim.run_timeblocks():
...     print(f'{block=}')  # Replace this with code to analyze the current configuration
block=0
block=1
block=2
block=3
block=4
block=5
block=6
block=7

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 run_timeblocks().

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.

• compute_flags (dict, optional) – Dictionary specifying which quantities to compute (volume, energies, stresses, etc.). If None the defaults are processed, see get_default_compute_flags():.

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)

The compute_flags paramter can be used if there should be allocated memory for storing volume data (or other data).

>>> configuration = gp.Configuration(D=3, compute_flags={'Vol':True})

The default values can be seen with get_default_compute_flags():

>>> gp.get_default_compute_flags()
{'U': True, 'W': True, 'K': True, 'lapU': False, 'Fsq': False, 'stresses': False, 'Vol': False, 'Ptot': False}

atomic_scale(density: float) → None

Scale atom positions and simulation box to a given density.

Parameters:

density (float) – Number density of particles after scaling.

copy_to_device() → None

Copy all data from host to device memory (CPU to GPU).

copy_to_host() → None

Copy all data from device to host memory (GPU to CPU).

classmethod from_h5(h5file: File, group_name: str, reset_images: bool = False, compute_flags: bool = None, include_topology: bool = True) → Configuration

Read a configuration from an open HDF5 file identified by group-name

Parameters:

• h5file (HDF5 File) – open HDF5 file object, as returned by h5py.File()

• group_name (str) – string corresponding to a key in the h5py.File containing a saved gamdpy configuration

• reset_images (bool) – if True set the images to zero (default False)

• compute_flags (bool) – NOTE: still to be developed, should be possible to define compute flags from dictionary compute_flags defining what will be stored in the configuration (default None)

• include_topology (bool) – if True then read also the topology from the file (default True)

Returns:

configuration – a gamdpy configuration object

Return type:

Configuration

Example

>>> import gamdpy as gp
>>> output_file = h5py.File('examples/Data/LJ_r0.973_T0.70_toread.h5')
>>> conf = Configuration.from_h5(output_file, 'restarts/restart0000')
>>> print(conf.D, conf.N, conf['r'][0])     # Print number of dimensions D, number of particles N and position of first particle
3 2048 [-6.407801 -6.407801 -6.407801]

get_potential_energy() → float

Get total potential energy of the configuration

get_volume() → float

Get volume of simulation box associated with configuration

make_lattice(unit_cell: dict, cells: list, rho: float = None, ptype_unit_cell: list = 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. The simulation box is Orthorhombic.

Unit cells directories are available in gamdpy.unit_cells.

Parameters:

• unit_cell (dict) – Dictionary with fractional_coordinates for particles and lattice_constants

• cells (list[int]) – Number of unit cells in each direction

• rho (float) – Number density

• ptype_unit_cell (list) – Types of the particles in the unit cell

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: float) → None

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. The simulation box type is Orthorhombic and cubic.

Parameters:

• N (int) – Number of particles

• rho (float) – Number density of particles

Example

>>> import gamdpy as gp
>>> configuration = gp.Configuration(D=3)
>>> configuration.make_positions(N=1000, rho=1.2)

randomize_velocities(temperature: float, seed=None, ndofs=None) → None

Randomize velocities according to a given temperature.

Parameters:

• temperature (float) – Temperature to randomize velocities by. If <= 0, set all velocities to zero.

• seed (int, optional) – Seed for random number generator

• ndofs (int, optional) – Number of degrees of freedom

Raises:

• ValueError – If spatial dimention (D) is None

• ValueError – If any mass is zero. Set masses before randomizing velocities.

save(output: File, group_name: str, mode: str = 'w', update_ptype: bool = True, include_topology: bool = True, use_topology_link: bool = False, verbose: bool = True) → None

Write a configuration to a HDF5 file

Parameters:

• configuration (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

• include_topology (bool) – Boolean flag indicating whether the topology of the configuration should be included

• use_topology_link (bool) – Boolean flag indicating the topology should just be a soft link to that of the initial configuration (if present)

• verbose (bool) – Boolean flag indicating whether messages should be written

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)

set_kinetic_temperature(temperature: float, ndofs=None) → None

Rescale velocities so magnitude correspond to a given temperature.

Parameters:

• temperature (float) – Temperature after rescaling of velocities.

• ndofs (int, optional) – Degrees of freedom. If not provided, ndofs = D*(N-1)

Raises:

ValueError – If the current temperature is zero (velocities are zero).

Simulation boxes

An simulation box object is attached to an configuration object.

class gamdpy.Orthorhombic(D: int, lengths: list)

Standard rectangular simulation box class

Parameters:

• D (int) – Spatial dimension

• lengths (list of floats) – The box lengths in each spatial dimension.

Raises:

ValueError – If the length of lengths is not equal to D.

Example

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

copy_to_device() → None

Copy SimulationBox data from host to device.

copy_to_host() → None

Copy SimulationBox data from device to host.

get_apply_PBC()

Apply periodic boundary conditions (PBC).

get_dist_moved_exceeds_limit_function()

For use in neighbor list : Single-particle criterion for whether the neighbor list needs to be rebuilt.

get_dist_sq_dr_function() → callable

Generates function dist_sq_dr which computes displacement and distance for one neighbor.

get_dist_sq_function() → callable

Generates function dist_sq_function which computes distance squared for one neighbor.

get_lengths() → ndarray

Return the box lengths as a numpy array

Returns:

The box lengths in each spatial dimension.

Return type:

numpy.ndarray

get_loop_x_addition()

For use in linked list implementation for neighbor list when Lees-Edwards (shearing) boundary conditions apply. In non-shearing cases zero should be returned.

get_loop_x_shift_function()

For use in linked list implementation for neighbor list when Lees-Edwards (shearing) boundary conditions apply. In non-shearing cases the function should return zero.

get_name() → str

Return name of the SimulationBox.

get_volume() → float

Return the box volume, \(V = \prod_{i=1}^{D} L_i\)

Returns:

Volume of the box.

Return type:

float

get_volume_function()

Returns the function which calculates the volume of the simulation box.

scale(scale_factor: float) → None

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.

Parameters:

• D (int) – Spatial dimension

• lengths (list of floats) – Lengths of the box sides.

• box_shift (float, optional) – Shift of the box in the x-direction (the direction of shearing). The default is no shift.

• box_shift_image (int) – Image shift of the box shift. The default is no shift.

Raises:

ValueError – If D is less than 2.

Example

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

copy_to_device()

Copy SimulationBox data from host to device.

copy_to_host()

Copy SimulationBox data from device to host.

get_apply_PBC()

Apply periodic boundary conditions (PBC).

get_dist_moved_exceeds_limit_function()

For use in neighbor list : Single-particle criterion for whether the neighbor list needs to be rebuilt.

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_lengths()

Return the box lengths as a numpy array

Returns:

The box lengths in each spatial dimension.

Return type:

numpy.ndarray

get_loop_x_addition()

For use in linked list implementation for neighbor list when Lees-Edwards (shearing) boundary conditions apply. In non-shearing cases zero should be returned.

get_loop_x_shift_function()

For use in linked list implementation for neighbor list when Lees-Edwards (shearing) boundary conditions apply. In non-shearing cases the function should return zero.

get_name()

Return name of the SimulationBox.

get_volume()

Return the box volume, \(V = \prod_{i=1}^{D} L_i\)

Returns:

Volume of the box.

Return type:

float

get_volume_function()

Returns the function which calculates the volume of the simulation box.

scale(scale_factor: float) → None

Scale the box lengths by scale_factor

Integrators

One of the below integrators should be given as a parameter to the Simulation class.

Constant energy and volume

class gamdpy.NVE(dt: float)

Total energy conserving integrator.

Consider a conservative force field \(f\). In this integrator, Newton’s equation of motion

\[m\ddot x = f\]

are discretized using the Leap-Frog algorithm

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

where \(v=\dot x\). This algorithm conserves the total energy up to numerical accuracy of floating point operations.

Parameters:

dt (float) – Time step for discretization

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

Constant temperature and volume

class gamdpy.NVT(temperature, tau: float, dt: float, integrator_state: ndarray = None)

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.

• integrator_state (Optional, zero-dimensional numpy array) – The value of the internal integrator state variable, relevant when restarting a simulation from a restart file

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

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

NVT Langevin Leap-frog integrator.

Leap-Frog implementation of the algorithm described in Ref. [Grønbech2014]. This integrator is a stochastic thermostat that keeps the system at a constant temperature, via the Langevin equations of motion

\[m \ddot x = f - \alpha \dot x + \beta\]

where \(f\) is the force from a conservative field, \(m\) is the particle mass, \(\alpha\) is a friction coefficient, and \(\beta\) is uncorrelated Gauss distributed noise: \(\langle \beta(t)\rangle=0\) and \(\langle \beta(t)\beta(t')\rangle=2\alpha T\delta(t-t')\) where \(T\) is the target temperature. Temperature is in reduced units where \(k_B=1\). For choosing the \(\alpha\) parameters, it is instructive to note that a characteristic timescale is given by

\[\tau = m/\alpha.\]

Parameters:

• temperature (float or function) – Temperature of the thermostat, \(T\). If a function, it must take a single argument, time, and return a float.

• alpha (float) – Friction coefficient of the thermostat, \(\alpha\).

• dt (float) – a time step for the integration.

References

[Grønbech2014]

Niels Grønbech-Jensen, Natha Robert Hayre, and Oded Farago, “Application of the G-JF Discrete-Time Thermostat for Fast and Accurate Molecular Simulations”, Comput. Phys. Commun. 185, 524-527 (2014) https://doi.org/10.1016/j.cpc.2013.10.006 https://arxiv.org/pdf/1303.7011.pdf

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

class gamdpy.Brownian(temperature, tau: float, dt: float, seed=0)

Brownian dynamics.

Implementation of Brownain dynamics (overdamped Langevin) using the algorithm GJ-0 described in Ref. [Grønbech2019] (Eq. (67)). This stochastic integrator keeps the system at a constant temperature, via the equations of motion

\[\alpha \dot x = f + \beta\]

where \(f\) is the force from a conservative field, \(\alpha\) is a friction coefficient, and \(\beta\) is uncorrelated Gauss distributed noise: \(\langle \beta(t)\rangle=0\) and \(\langle \beta(t)\beta(t')\rangle=2\alpha T\delta(t-t')\) where \(T\) is the temperature of a solvent. In this implimentatio, the friction coefficient \(\alpha\) is given by

\[\alpha = m/\tau\]

where \(m\) is the mass of a given particle, and \(\tau\) is a characteristic time. Thus, the particle mass, \(m\), can be used to set individual particle couplings with the solvent.

The equation of motion is discretized using a time step \(dt\) by

\[x(t+dt) = x(t) + \frac{1}{\alpha}f(t)dt+\frac{1}{2\alpha}[\beta(t)+\beta(t+dt)]\]

The “current” velocity stored is defined as

\[v(t) = [x(t-dt) - x(t)]/dt\]

Parameters:

• temperature (float or function) – Temperature of the implicit solvent. If a function, it must take a single argument, time, and return a float.

• tau (float) – Collision time of implicit solvent.

• dt (float) – a time step for the integration.

• seed (int, optional) – Seed for the pseudo random \(\beta\)’s.

References

[Grønbech2019]

Niels Grønbech-Jensen, “Complete set of stochastic Verlet-type thermostats for correct Langevin simulations”, Mol. Phys. 118, e1662506 (2019) https://doi.org/10.1080/00268976.2019.1662506

Examples

>>> import gamdpy as gp
>>> integrator = gp.Brownian(temperature=1.0, tau=0.1, dt=0.005)

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

Constant temperature and pressure

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

Constant NPT integrator for atomic systems.

Integrator that keeps the number of particles (\(N\)), pressure (\(P\)), and temperature (\(T\)) constant, using a Leap-Frog implementation of the Nosé–Hoover thermostat by Martyna–Tuckerman–Klein presented in Ref. [Martyna1996]. (Note that the thermostat and barostat states in this implementation are defined as \(p_\xi/Q\) and \(p_\epsilon/W\), see Eq. 2.9 in the reference.).

The thermal and barometric coupling parameters are defined via two time-scales.

Parameters:

• temperature (float or function) – Target temperature

• tau (float) – Thermostat relaxation time

• pressure (float or function) – Target pressure

• tau_p (float) – Barostat relaxation time

• dt (float) – a time step for the integration.

Raises:

• TypeError – If the simulation box is not Orthorhombic.

• ValueError – If the spatial dimension of the simulation box is not \(D=3\).

References

[Martyna1996]

Glenn J. Martyna, Mark E. Tuckerman, Douglas J. Tobias, and Michael L. Klein, “Explicit Reversible Integrators for Extended Systems Dynamics”, Mol. Phys. 87, 1117–57 (1996) https://doi.org/10.1080/00268979600100761

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

class gamdpy.NPT_Langevin(temperature, pressure, alpha: float, alpha_barostat: float, mass_barostat: float, dt: float, seed=0, integrator_state: ndarray = None)

Constant NPT Langevin integrator with isotropic volume fluctuations.

This integrator is a Leap-Frog implementation of the GJ-F Langevin equations of motion present in Ref. [Grønbech2014b]. It is intended for an orthorhombic simulation box in any spatial dimention.

Let \(f\) be the conservative force field, \(\alpha\) is a friction parameter, and \(\beta\) be uncorrelated Gauss distributed noise: \(\langle \beta(t)\rangle=0\) and \(\langle \beta(t)\beta(t')\rangle=2\alpha k_B T\delta(t-t')\) where \(T\) is the target temperature. The Langevin equation of motion for a given particle is

\[m\ddot x = f - \alpha \dot x + \beta\]

For choosing the \(\alpha\) parameters, it is instructive to note that a characteristic timescale is given by

\[\tau_T = m/\alpha\]

The (isotropic) equations of motions for the volume of the simulation box (\(V\)) can similarly be written as

\[Q\ddot V = f_P - \alpha_V \dot V + \beta_V\]

Here, \(Q\) is an inertial coefficient (a barostat mass), \(\alpha_V\) is another friction parameter, \(\beta_V\) is uncorrelated Gauss distributed noise, and \(f_P\) is a “barostat force”. Specifically, if \(W\) is the (instantaneous) virial, then

\[f_P = W + \frac{Nk_B T}{V} - P\]

where \(P\) is the target pressure. To set the barstat parameters \(Q\) and \(\alpha_V\) it can be instructive to note that equations of motions for the volume resemble that of a damped harmonic oscillator. If \(K=-VdP/dV\) is the bulk modulus of the system, then the “spring constant” of the oscillator is \(k=K/V\), then the natural frequency is \(\omega_0=\sqrt{k/Q}\) and damping ratio is \(\zeta=\alpha_V/2\sqrt{Qk}\) (\(\zeta<1\) for underdamped oscillation). If a characteristic timescale is defined as \(\tau_V\equiv 1/\omega_0\), then

\[\alpha_V = 2 K \tau_v/V\]

\[Q = K (\zeta \tau_v)^2 / V\]

Parameters:

• temperature (float or function) – Target temperature, \(T\)

• pressure (float or function) – Target pressure, \(P\)

• alpha (float) – Friction coefficient of the thermostat, \(\alpha\)

• alpha_barostat (float) – Friction coefficient of the barostat, \(\alpha_V\)

• mass_barostat (float) – Inertial coefficient (barostat mass), \(Q\)

• dt (float) – a Time step for the Leap-Frog discretization

• volume_velocity (float) – Initial velocity of volume

• seed (int) – seed for the (pseudo) random noise

Raises:

TypeError – If the simulation box is not Orthorhombic.

References

[Grønbech2014b]

Niels Grønbech-Jensen and Oded Farago, “Constant pressure and temperature discrete-time Langevin molecular dynamics”, J. Chem. Phys. 141, 194108 (2014) https://doi.org/10.1063/1.4901303

Examples

Example of setting parameters for the NPT langevin barostat (see above) in “Lennard-Jones like” reduced units.

>>> import gamdpy as gp
>>> tau_T = 2.0
>>> tau_V = 8.0
>>> zeta = 0.2
>>> K = 100   # Replace with an estimate of the bulk modulus
>>> V = 1000  # Replace with an estimate of the average volume
>>> integrator = gp.NPT_Langevin(
...    temperature=2.0,
...    pressure=1.0,
...    alpha=1/tau_T,
...    alpha_barostat=2*K/tau_V/V,
...    mass_barostat=K*(zeta*tau_V)**2/V,
...    dt=0.004)

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

Other integrators

class gamdpy.SLLOD(shear_rate, dt)

The SLLOD integrator

Shear an atomic system in the xy-plane using the SLLOD equations [Edberg1986] implimented with the operator splitting algorithm, see [Pan2005]. This integrator needs a LeesEdwards simulation box.

Parameters:

• shear_rate (float) – The shear rate of the system.

• dt (float) – The time step of the simulation.

Raises:

TypeError – If the simulation box is not LeesEdwards simulation box.

References

[Edberg1986]

Roger Edberg, Denis J. Evans and G. P. Morriss

“Constrained molecular dynamics: Simulations of liquid alkanes with a new algorithm” J. Chem. Phys. 84, 6933–6939 (1986) https://doi.org/10.1063/1.450613

[Pan2005]

Guoai Pan, James F. Ely, Clare McCabe and Dennis J. Isbister

“Operator splitting algorithm for isokinetic SLLOD molecular dynamics” J. Chem. Phys. 122, 094114 (2005) https://doi.org/10.1063/1.1858861

Examples

An example of how to set up a Lees Edwards simulation box and a SLLOD integrator in a Lennard-Jones like system.

>>> import gamdpy as gp
>>> configuration = gp.Configuration(D=3, compute_flags={'stresses': True})
>>> configuration.make_lattice(gp.unit_cells.FCC, cells=[8, 8, 8], rho=0.973)
>>> configuration['m'] = 1.0
>>> configuration.randomize_velocities(temperature=0.7)
>>> configuration.simbox = gp.LeesEdwards(configuration.D, configuration.simbox.get_lengths())
>>> configuration.set_kinetic_temperature(temperature=0.7, ndofs=configuration.N*3-4)  # Note that we need to subtract 4
>>> integrator = gp.SLLOD(shear_rate=0.02, dt=0.01)

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

class gamdpy.GradientDescent(dt: float)

Gradient descent algorithm, minimizing the potential energy.

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

Parameters:

dt (float) – Time step for discretization / Learning rate

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

class gamdpy.ActiveOUP(DT: float, DA: float, mu: float, tau: float, dt: float, seed=0)

Active Ornstein-Uhlenbeck Particle (AOUP).

Implementation of overdamped active Ornstein-Uhlenbeck dynamics. Active Ornstein-Uhlenbeck Particle, is an active matter model in which overdamped particles are subject to a colored Ornstein-Uhlenbeck noise, see Ref. [Fodor2016].

The equations of motion:

\[\dot{x_i} = -\mu f_i + \eta_i + \xi_i\]

where \(\mu\) is a mobility, and \(f\) is an interaction potential. \(\xi\) is an uncorrelated gauss distributed thermal noise:

\[\langle \xi_{i_\alpha}(t)\xi_{j\beta}(t')\rangle = 2D_T \delta_{ij} \delta_{\alpha \beta} \delta(t-t')\]

here \(D_T\) is a diffusion coefficient, and Greek letters correspond to spatial components. \(\eta_i\) are Ornstein-Uhlenbeck processes, solution of

\[\dot{\mathbf{\eta}}_i = -\frac{\mathbf{\eta}_i}{\tau} + \frac{\sqrt{D_A}}{\tau} \mathbf{\zeta}_i\]

with \(\langle \zeta_{i\alpha}(t) \zeta_{j\beta}(t')\rangle = 2\delta_{ij} \delta_{\alpha \beta} \delta(t-t')\)

The autocorrelation of \(\eta\) then is

\[\langle \eta_{i\alpha}(t)\eta_{j\beta}(t')\rangle= \delta_{ij} \delta_{\alpha \beta} \frac{D_A}{\tau} e^{\frac{-|t-t'|}{\tau}}\]

so \(D_A\) controls the amplitude of the noise and \(\tau\) its persistence time.

The equations are discretized using the Euler-Maruyama method [Higham2001] with a timestep \(dt\):

\[\mathbf{r}_i(t+dt) = \mathbf{r}_i(t)+(\mu\mathbf{F}_i + \mathbf{\eta}_i)dt + \sqrt{2 D_T}\sqrt{dt} N(0,1)\]

\[\mathbf{\eta}_i(t+dt) = \mathbf{\eta}_i(t)-\frac{\mathbf{\eta}_i}{\tau}dt + \frac{\sqrt{2 D_A}}{\tau} N(0,1)\]

where \(N(0,1)\) is a normal distributed pseudo random number.

Parameters:

• DT (float) – Thermal diffusion coefficient

• DA (float) – “Active” diffusion coefficient

• mu (float) – mobility

• tau (float) – persistence time of the colored noise

• dt (float) – timestep

• seed (int, optional) – seed for pseudo random number generator

References

[Fodor2016]

Etienne Fodor et al. “How Far from Equilibrium Is Active Matter?” Phys. Rev. Lett. 117 (3 July 2016), p. 038103 https://doi.org/10.1103/PhysRevLett.117.038103

[Higham2001]

Higham, D. J. (2001). “An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations.” SIAM Review, 43(3), 525–546. https://doi.org/10.1137/S0036144500378302

Examples

>>> import gamdpy as gp
>>> integrator = gp.integrators.ActiveOUP(DT=0.25, DA=0.5, mu=1.0, tau=0.75, dt=0.001, seed=0)

save_internal_state(output: File, group_name: str)

Write the internal state of the integrator as an attribute called ‘integrator_state’ to the specified group in the specified HDF file

Interactions

Interactions are passed in a list to the Simulation class. This will typically include a pair potential and fix interactions like gravity and walls.

Pair potential

Evaluating pair potentials is typically computationally expensive. For efficiency, instances of the class below only compute forces between two particles when their separation is less than a specified cutoff distance.

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

Pairwise interaction potential for a system of particles.

Parameters:

• pairpotential_function (callable) –

  A JIT compiled function f(r, params) that takes a separation distance r (float) and a list of parameters, and returns a triplet of floats:

  • Potential energy, \(u(r)\)

  • Force multiplier, \(-u'(r)/r\)

  • Second derivative of potential energy, \(u''(r)\)

• params (list of floats or nested list of floats) – Interaction parameters for the pair potential function. Use a nested list for multiple types of particles. The last element of each list is the cutoff distance of the pair potential.

• max_num_nbs (int) – Maximum number of neighbors per particle to allocate in the neighbor list.

• exclusions (array_like) – List of particle indices to exclude from interactions for each particle.

Example

The standard Lennard-Jones potential shifted and truncated at 2.5:

>>> pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon)
>>> sig, eps, cut = 1.0, 1.0, 2.5
>>> pair_pot = gp.PairPotential(pair_func, params=[sig, eps, cut], max_num_nbs=1000)
>>> interactions = [pair_pot, ]  # Passed to a Simulation instance

merge_exclusions(exclusion_list, max_exclusions=None)

Take a list of exclusion arrays are return a merged array

Pair potential functions

gamdpy.LJ_12_6(dist, params)

The 12-6 Lennard-Jones potential

\[u(r) = A_{12} r^{-12} + A_6 r^{-6}\]

See gamdpy.apply_shifted_potential_cutoff() for a usage example of a shifted potential cutoff.

Parameters:

• dist (float) – Distance between particles

• params (array-like) – \(A_{12}\), \(A_{6}\)

Returns:

• u (float) – Potential energy, \(u(r)\)

• s (float) – Force multiplier, \(-u'(r)/r\)

• umm (float) – Second derivative of the potential energy, \(u''(r)\)

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.yukawa(dist, params)

The Yukawa potential (simple screened Coulomb potential)

\[u(r) = \varepsilon \frac{\sigma}{r} \exp(-r/\sigma)\]

Parameters:

• dist (float) – Distance between particles

• params (array-like) – σ, ε

gamdpy.gaussian_core_model(dist, params)

Gaussian-core model (GCM) [Stillinger1976]

\[u(r) = \varepsilon \exp(-r^2/\sigma^2)\]

Parameters:

• dist (float) – Distance between particles

• params (array-like) – [sigma, epsilon]

References

[Stillinger1976]

Frank H. Stillinger, “Phase transitions in the Gaussian core system”, J. Chem. Phys. 65, 3968–3974 (1976), https://doi.org/10.1063/1.432891

gamdpy.exponential_repulsion(dist, params)

The exponential repulsion pair potential (EXP)

\[u(r) = \varepsilon \exp(-r/\sigma)\]

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₅, σ, ε

gamdpy.universal_zbl_potential(dist, params)

The Ziegler-Biersack-Littmark (ZBL) universal potential for high-energy pair collisions of atoms.

The ZBL potential is a sum of four screened Coulomb interactions [Ziegler1985] :

\[u(r) = \varepsilon \sum^4_{k=1} \frac{\sigma}{r} c_k\exp(-b_k r/\sigma)\]

where the c’s are [0.18175, 0.50986, 0.28022, 0.02817] and the b’s are [3.19980, 0.94229, 0.40290, 0.20162]. These are the same values as used in the LAMMPS implementation.

Consider to atoms with atomic numbers \(Z_i\) and \(Z_j\). Then the screening length is

\[\sigma = \frac{0.46850 \text{Å}}{Z_i^{0.23}+Z_j^{0.23}}\]

and the energy parameter is

\[\varepsilon = \frac{Z_i Z_j e^2}{4 \pi \epsilon_0 \sigma}\]

As an example, for copper (\(Z_i=Z_j=29\)) the screening length is \(\sigma=0.1080\) Å, the energy parameter is \(\varepsilon=1.122\times10^5\) eV, \(\varepsilon/k_B=1.301\times10^9\) K.

Parameters:

• dist (float) – Distance between particles

• params (array-like) – \(\sigma\), \(\varepsilon\)

References

[Ziegler1985]

JF Ziegler, JP Biersack, U Littmark “The Stopping and Range of Ions in Solids”, Pergamon, 1985

Generators

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

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

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

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_m(r) = u(r) - u(r_{c}) \quad (r<r_c)\]

and \(u_m(r)=0\) for \(r>r_c\) where \(r_c\) is the cutoff distance. Calls the original potential function twice, avoiding changes to params. The last entry the parameter array (params) is the cutoff: […, r_c].

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 potentia cutoff is applied to the original function.

Return type:

callable

Example

Example demonstrating how to set up the Lennard-Jones 12-6 potential, LJ_12_6(), truncated and shifted to zero at the cutoff distance of 2.5.

>>> import gamdpy as gp
>>> pair_func = gp.apply_shifted_potential_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)
>>> interactions = [pair_pot, ]  # List of interactions passed to an instance of the Simulation class

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_m(r) = u(r) - u(r_{c}) + s(r_{c})(r - r_{c}) \quad (r<r_c)\]

and \(u_m(r)=0\) for \(r>r_c\), where \(r_c\) is the cutoff distance, and \(s(r) = -u'(r)/r\). The last entry the parameter array (params) is the cutoff: […, r_c].

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

gamdpy.apply_cubic_spline_cutoff(pair_potential)

Apply cubic spline cutoff to a pair-potential function

Actually the cubic spline is applied to the pair force as done by LAMMPS (see https://docs.lammps.org/pair_lj_smooth.html). The potential energy is given by a quartic function in the transition region. Specifically, if \(u(r)\) is the original potential, then the (modified) truncated potential is

\[u_m(r) = u(r) + K \quad (r < r_i),\]

\[u_m(r) = C_0 + C_1 [r-r_i] + \frac{C_2}{2} [r-r_i]^2 + \frac{C_3}{3} [r-r_i]^3 + \frac{C_4}{4} [r-r_i]^4 \quad (r_i \le r \le r_c),\]

and

\[u_m(r) = 0 \quad (r > r_c),\]

The parameters (the \(C\)’s and the \(K\)) for a smooth truncation are computed automatically: At the inner cutoff radius, the force and its first derivative match the unsmoothed potential (\(u(r)\)). At the outer cutoff radius, both the force and the first derivative are zero. See SM of Ref. [Leoni2025] for more details.

The last two values of the parameter array (params) are the inner and outer cutoffs: […, r_i, r_c].

Parameters:

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

Returns:

potential – a function where a cubic spline cutoff is applied to original function

Return type:

callable

References

[Leoni2025]

Fabio Leoni, John Russo, Francesco Sciortino, and Taiki Yanagishima. “Generating Ultrastable Glasses by Homogenizing the Local Virial Stress”. Phys. Rev. Lett. 134, 128201 (2025) https://doi.org/10.1103/PhysRevLett.134.128201

gamdpy.apply_gromacs_cutoff(pair_potential: callable) → callable

Apply a smooth cutoff to a pair-potential function using GROMACS style.

The potential is smoothly truncated at \(r_c\) by applying a switching function. Specifically, if \(u(r)\) is the original potential, then the (modified) truncated potential is

\[u_m(r) = u(r) + S(r)\]

Outside the cut-off, the switching function is \(S(r) = -u(r)\) (\(r > r_c\)) truncating the potential at \(r_c\) \(u_m(r) = 0 \quad (r > r_c)\). In the switching region

\[S(r) = \frac{A}{3}[r-r_1]^3 + \frac{B}{4}[r-r_1]^4 + C \quad (r_1 < r < r_c)\]

At \(r<r_1\), the potential is shifted by a constant

\[S(r) = C \quad (r < r_1)\]

The parameters ensuring smooth truncation are

\[A = \frac{- 3 u'(r_c) + [r_c - r_1] u''(r_c) }{[ r_c - r_1 ]^2}\]

\[B = [ 2 u'(r_c) - [r_c - r_1] u''(r_c) ]/[ r_c - r_1 ]^3\]

\[C = - u(r_c) + \frac{1}{2} [r_c - r_1] u'(r_c) - \frac{1}{12} [r_c - r_1]^2 u''(r_c)\]

The parameters are computed automatically. The last two values of the parameter array (params) are the inner and outer cutoffs: […, r_1, r_c].

For more, see GROMACS force-switch function (3rd degree polynomial) in manual at https://www.gromacs.org/, or the LAMMPS documentation, https://docs.lammps.org/pair_gromacs.html.

Parameters:

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

Returns:

potential – a function where a smooth cutoff is applied to the original function.

Return type:

callable

Example

>>> import gamdpy as gp
>>> params = sig, eps, cut_inner, cut_outer = 1.0, 1.0, 3.0, 4.0
>>> pair_func = gp.apply_gromacs_cutoff(gp.LJ_12_6)

Fixed interactions

Fixed interactions (covalent bonds, angles, gravitational forces, walls, or tethers to an anchor points) are evaluated at every time step. This is contrast to non-bonded pair interactions, which are only computed for particle pairs within the non-bonded cutoff distance.

Classes

class gamdpy.Bonds(bond_potential, bond_indices, bond_parameters)

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.

• bond_indices (list) – A list of lists, each containing the indices of the two particles involved in a bond and the bond

• bond_parameters (list) – A list of parameters for each bond type. Each entry is a list of parameters for a specific bond type.

See also

gamdpy.harmonic_bond_function

Harmonic bond potential

class gamdpy.Angles(angle_potential, angle_indices, angle_parameters)

Fixed angle interactions between particle triplet

Parameters:

• angle_potential (function) – A angle potential function, see gamdpy.harmonic_angle_function() for an example.

• indices (list) – A list of indices, each containing the indices of the particles involved in a bond and the bond type.

• potential_params (list) – A list of parameters for each angle type.

See also

gamdpy.harmonic_angle_function

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.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)

class gamdpy.Planar(potential: Callable, params: list[list[float]], indices: list[list[int]], normal_vectors: list[list[float]], points: list[list[float]])

Planar interactions such as smooth walls, gravity, or an electric field.

Consider a plane with the normal vector \({\bf n}\) going though the point \({\bf p}\). For a given particle, let \({\bf r}\) be the distance to the nearest point in the plane. Then the planar force is

\[{\bf F} = s(r) {\bf n}\]

Where \(s(r)=-u'(r)/r\) is the force multiplier of a given potential function.

Note: The planer interaction is considered an “external force”, and will not contribute particles scalar energies, virials etc.

Parameters:

• potential (Callable) – Potential function for planer interactions See :func:gamdpy.potential_functions.harmonic_repulsion for an example.

• params (list[list[float]]) – A list of parameters for each plane type. Each entry is a list of parameters for the potential function (see above).

• indices (list[list[int]]) – A list of lists, each containing the indices of a particle involved, and the planer interactions of relevance.

• normal_vectors (list[list[float]]) – A list of lists, each containing a normal vector for a given plane.

• points (list[list[float]]) – A list of lists, each containing a point on a given plane.

Example

Planar interactions can be used to add smooth walls and gravity to a simulation (be aware of periodic boundary conditions).

>>> import gamdpy as gp
>>>
>>> # Replace the below with your own simulation
>>> sim = gp.get_default_sim()
>>> box_length = sim.configuration.simbox.get_lengths()[1]  # Box-length in y-direction
>>> N = sim.configuration.N
>>>
>>> # Two smooth walls
>>> wall_distance = box_length/2
>>> walls = gp.interactions.Planar(
...     potential=gp.harmonic_repulsion,
...     params=[[100.0, 1.0], [100.0, 1.0]],
...     indices=[[n, 0] for n in range(N)] + [[n, 1] for n in range(N)],  # All particles feel both walls
...     normal_vectors=[[0.0, 1.0, 0.0], [0.0, -1.0, 0.0]],
...     points=[[0.0, -wall_distance/2, 0], [0.0, wall_distance/2, 0]]
... )
>>>
>>> # Gravity
>>> mg = 0.0005
>>> potential_gravity = gp.make_IPL_n(-1)
>>> gravity = gp.interactions.Planar(
...     potential=potential_gravity,
...     params= [[mg, 10*wall_distance]],
...     indices= [[n, 0] for n in range(N)],   # All particles feel the gravity
...     normal_vectors= [[0,1,0], ],
...     points= [[0, -wall_distance/2.0, 0] ]  # Defining 0 for potential energy on the lower wall
... )
>>> interactions = [..., walls, gravity]

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.

Bond functions

A bond potential function is needed for the Bonds class.

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

Angle functions

An angle potential function is needed for the Angles class.

gamdpy.harmonic_angle_function(theta: float, params: ndarray) → tuple

Harmonic angle potential,

With a regularization parameter SMALL = 1.e-6 hard-coded in. To change the value of SMALL call make_harmonic_angle_function(SMALL) .. math:

u(\theta) = \frac{k}{2} (\theta - \theta_0)^2

Parameters:

• theta (float) – Angle (radians) defined by three neighboring particles in a molecule, more precisely: the angle subtended by atoms 0 and 2 at 1

• params (array-like) – \(\theta_0\), the angle of minimum energy, \(k_{spring}\), the spring constant.

Returns:

• u (float) – Potential energy

• d_u_cos_theta_neg (float) – Negative derivative of the potential energy with respect to cos(theta). This involves dividing by sin(theta)+SMALL

See also

gamdpy.Angles

gamdpy.cos_angle_function(angle, params)

Cosine angle potential,

\[u(\theta) = \frac{k}{2} (\cos(\theta) - \cos(\theta_0))^2\]

Parameters:

• theta (float) – Angle (radians).

• params (array-like) – [theta_0, k] Angle of minimum energy (zero force)) and spring constant.

See also

gamdpy.harmonic_angle_function

Runtime Actions

A list of runtime actions are passed as an argument to the Simulation class.

class gamdpy.TrajectorySaver(scheduler=<gamdpy.runtime_actions.time_scheduler.Log2 object>, include_simbox=False, compression='gzip', compression_opts=4, saving_list=['positions', 'images'], update_ptype=False, update_topology=False, verbose=False)

Runtime action for saving configurations during timeblock.

Parameters:

• scheduler (Scheduler) – The scheduler defining when to save

• include_simbox (bool, optional) – Boolean deciding if to include simbox informations in output Default False

• compression (str, optional) – String selecting the type of compression used. Default “gzip”

• compression_opts (int, optional) – Relevant if “gzip” compression is used. Select compression option for gzip

• saving_list (list, optional) – This list is used to select which information to save in the trajectory. Default [‘positions’, ‘images’]. Can be used to save only velocities by setting saving_list = [‘velocities’]. Possible list entries are ‘positions’, ‘images’, ‘velocities’, ‘forces’.

• verbose (bool, optional) – Default False

Raises:

ValueError – If saving_list contains a name which is not recongnized.

class gamdpy.ScalarSaver(steps_between_output: int = 0, 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.

columns() → list[str]

Get a list of available data columns saved by ScalarSaver in a .h5 file

Parameters:

h5py.File (h5file) – HDF5 file object from which data is read

Returns:

list(str)

Return type:

A list of available data column keys

extract(columns: list[str], per_particle: bool = True, first_block: int = 0, last_block: int = None, subsample: int = 1, function: callable = None) → list

Get a tuple of lists of time series of available data columns saved in a .h5 file

Parameters:

• h5file (h5py.File) – HDF5 file object from which data will be read

• columns (list[str]) – List of keys for data columns to be extracted, eg. [‘U’, ‘K’,]

• per_particle (bool) – Bolean flag determining whether date should be returned divided by number of particles (default True)

• first_block (int) – First timeblock to include in returned data (default 0)

• last_block (int or None) – last timeblock to include in returned data (default None, i.e. include last available timeblock)

• subsample (int) – If ‘2’ return every second entry in saved time-series, etc (default 1)

• function (callable) – function applied to each time series data before returning, e.g.: np.mean (default None)

Returns:

list

Return type:

A list numpy arrays, one for each column asked for - or the results of applying ‘function’ to these

extract_as_dict(**kwargs)

Get a dictionary of time series of available data columns saved in a .h5 file

Same as ScalarSaver.extract(), but returns a dictionary. See ScalarSaver.extract() for more details on the arguments.

If ‘columns’ is not specified in kwargs, all available columns are returned.

get_times(first_block=0, last_block=None, reset_time=True, subsample=1)

Get a numpy array with times associated with data columns saved by ScalarSaver in a .h5 file

Parameters:

• h5file (h5py.File) – HDF5 file object from which data will be read

• first_block (int) – First timeblock to include in returned data (default 0)

• last_block (int ot None) – last timeblock to include in returned data (default None, i.e. include last available timeblock)

• reset_time (bool) – Bolean flag determining whether the returned times should start from zero (default True)

• subsample (int) – If ‘2’ return every second entry in saved time-series, etc (1)

Returns:

• np.array (A numpy arrays with simulation times)

• Seed also

• ———

• :function:`ScalarSaver.extract`

info() → str

Get information about what data was saved by ScalarSaver in a .h5 file

Parameters:

h5py.File (h5file) – HDF5 file object from which data is read

Returns:

str

Return type:

A string containing information about the available scalars, and the associated shape.

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:

‘distances’ - numpy array with distances to the middle of the bins ‘rdf_per_frame’ - numpy array [bin_index, typeA, typeB, frame] ‘rdf’ - numpy array [bin_index, typeA, typeB], i.e. averaged over frames ‘ptype’ - numpy array with particle types

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/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/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, mode='w')

This method saves self.h5 to disk. It can be used to save sim.output to file if TrajectoryIO class is initialized without arguments and then using self.h5 = sim.output.

LC: By default each runtime action has a compression setting. Would be nice if this function can change that when saving.

IO functions

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_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_to_lammps(configuration, timestep=0, unit_rescale=None) → str

Convert a configuration to a string formatted as LAMMPS dump file

Parameters:

• configuration (gamdpy.Configuration) – a gamdpy configuration object

• timestep (int) – time at which the configuration is saved

• unit_rescale (dictionary) – contains rescale factors for positions/box, velocities and forces for converting units.

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.thermodynamics_NpT(N: int, dof: int, U: Iterable[float], W: Iterable[float], K: Iterable[float], Vol: Iterable[float], k_B=1.0, T_ext: float = None, p_ext: float = None, per_particle=True) → dict

Calculate thermodynamic response functions from equilibrium fluctuations in an isotropic NpT simulation.

The implementation analyses time series of thermodynamic observables from a constant \(NpT\) (isothermal–isobaric) ensemble, specifically potential energy \(U\), configurational virial \(W= \left. \frac{\partial U\!\left(\lambda\mathbf{r}^N\right)}{\partial \lambda} \right|_{\lambda=1}\), kinetic energy \(K\), and volume \(V\), and computes thermodynamic response functions. The implementation assumes:

• An isotropic \(NpT\) ensemble at equilibrium.

• If per_particle=True (default), the inputs U, W, K, Vol are per-particle values and if per_particle=False, the inputs are treated as values for the entire system. N is used in the conversion.

• The external temperature, \(T\), defaults to \(\langle T_\mathrm{inst}\rangle\) if not provided, where the instantaneous temperature is computed from equipartition, \(T_\mathrm{inst} = \frac{2K}{k_B\, N_d}\) and \(\mathrm{N_d}\) is the number of degrees of freedom (dof).

• The external pressure, \(p\), defaults to \(\langle p_\mathrm{inst}\rangle\) if not provided, where the instantaneous pressure follows from, \(p_\mathrm{inst} = \frac{Nk_B T_\mathrm{inst} + W}{V}\).

The isobaric heat capacity (per particle), \(c_p=\frac{1}{N}\left(\frac{\partial H}{\partial T}\right)_p\), the isothermal compressibility, \(\kappa_T=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_T\), and the isobaric expansion coefficient, \(\alpha_p = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_p\) are computed using fluctuation-disipation formulas of the NpT ensemble. If \(H = U + K + p V\) is the instantaneous enthalpy, then

\[\begin{split}c_p &= \frac{\mathrm{Var}(H)}{N\, k_B\, T^2}, \\ \kappa_T &= \frac{\mathrm{Var}(V)}{\langle V\rangle\, k_B\, T}, \\ \alpha_p &= \frac{\mathrm{Cov}(V,H)}{\langle V\rangle\, k_B\, T^2}.\end{split}\]

From these, other thermodynamic identities are computed (see Returns below).

Parameters:

• N (int) – Number of particles, \(N\).

• dof (int) – Total number of degrees of freedom of the system, \(N_d\).

• U (Iterable[float]) – Time series of potential energy, \(U\).

• W (Iterable[float]) – Time series of configurational virial, \(W\).

• K (Iterable[float]) – Time series of kinetic energy, \(K\).

• Vol (Iterable[float]) – Time series of volume, \(V\).

• k_B (float, optional) – Boltzmann constant of the unit system (default 1.0), \(k_B\).

• T_ext (float, optional) – External (bath) temperature, \(T\). If None, then \(\langle T_\mathrm{inst}\rangle\) is used.

• p_ext (float, optional) – External (bath) pressure, \(p\). If None, then \(\langle p_\mathrm{inst}\rangle\) is used.

• per_particle (bool, optional) – If True (default), the U, W, K, Vol inputs are interpreted as per-particle values. If False, they are treated as extensive.

Returns:

A dictionary with scalar properties and response functions reported as intensive.

• density: \(\rho = N \langle 1/V \rangle\).

• specific_volume: \(v = \langle V \rangle/N\).

• jensen_bias_density_volume: \(\rho v - 1\)

• potential_energy: \(\langle U \rangle / N\).

• kinetic_energy: \(\langle K \rangle / N\).

• internal_energy: \(\langle E \rangle / N\) where \(E=U+K\).

• kinetic_temperature: \(\langle T_\mathrm{inst} \rangle\), with \(T_\mathrm{inst} = 2K/k_B N_d\).

• external_temperature: Provided external \(T\) (if given).

• pressure: \(\langle p_\mathrm{inst} \rangle\), with \(p_\mathrm{inst} = (N k_B T_\mathrm{inst} + W)/V\).

• external_pressure: Provided external \(p\) (if given).

• compressibility_factor: \(Z = p/\rho\, k_B T\).

• enthalpy: \(\langle H \rangle / N\), with \(H = U + K + p V\).

• isobaric_heat_capacity: \(c_p = \mathrm{Var}(H)/N k_B T^2\).

• isothermal_compressibility: \(\kappa_T = \mathrm{Var}(V) / \langle V \rangle k_B T\).

• isothermal_bulk_modulus: \(K_T = 1/\kappa_T\).

• isobaric_expansion_coefficient: \(\alpha_p = \mathrm{Cov}(V,H) / \langle V \rangle k_B T^2\).

• isochoric_heat_capacity: \(c_V = \frac{1}{N}\left(\frac{\partial E}{\partial T}\right)_V = c_p - \dfrac{T\, \alpha_p^2}{\rho\, \kappa_T}\).

• isochoric_heat_capacity_excess: \(c_V^{\mathrm{ex}} = \frac{1}{N}\left(\frac{\partial U}{\partial T}\right)_V = c_V - c_V^{\mathrm{id}}\), with \(c_V^{\mathrm{id}} = \frac{1}{N}\left(\frac{\partial K}{\partial T}\right)_V = \dfrac{N_d}{2N} k_B\).

• adiabatic_index: \(\gamma = \dfrac{c_p}{c_V}\).

• adiabatic_compressibility: \(\kappa_S = -\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_S = \kappa_T \dfrac{c_V}{c_p}\).

• adiabatic_bulk_modulus: \(K_S = 1/\kappa_S\).

• isochoric_pressure_coefficient: \(\beta_v = \left(\dfrac{\partial p}{\partial T}\right)_V = \dfrac{\alpha_p}{\kappa_T}\).

• isochoric_pressure_coefficient_excess: \(\beta_v^{\mathrm{ex}} = \left(\dfrac{\partial W}{\partial T}\right)_V / V =\beta_v-\rho k_B\).

• adiabatic_pressure_coefficient: \(\beta_s = \left(\dfrac{\partial p}{\partial T}\right)_S = \dfrac{\rho\, c_p}{T\, \alpha_p}\).

• adiabatic_expansion_coefficient: \(\alpha_s = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_s = \alpha_p - \kappa_T\, \beta_s\).

• thermodynamic_gruneisen_parameter: \(\gamma_G = V\left(\frac{\partial p}{\partial E}\right)_V = \left(\frac{\partial \ln T}{\partial \ln \rho}\right)_S = \dfrac{\alpha_p}{\rho \, \kappa_T \, c_V}\).

• configurational_adiabatic_scaling_exponent: \(\gamma = \left(\dfrac{\partial \ln T}{\partial \ln \rho}\right)_{S_\textrm{ex}} = \beta_v^{\mathrm{ex}} / \rho c_V^{\mathrm{ex}}\).`

• joule_thomson_coefficient: \(\mu_{JT} = \left(\dfrac{\partial T}{\partial p}\right)_H = \dfrac{T\, \alpha_p - 1}{\rho\, c_p}\).

Return type:

dict

gamdpy.tools.thermodynamics_NVT(N: int, dof: int, V: float, U: Iterable[float], W: Iterable[float], K: Iterable[float], k_B=1.0, T_ext: float = None, per_particle=True) → dict

Calculate thermodynamic response functions from equilibrium fluctuations in an isotropic NVT simulation.

The implementation analyzes time series of thermodynamic observables from a constant \(NVT\) (canonical) ensemble, specifically potential energy \(U\), configurational virial \(W= \left. \frac{\partial U\!\left(\lambda\mathbf{r}^N\right)}{\partial \lambda} \right|_{\lambda=1}\), kinetic energy \(K\), and computes thermodynamic response functions. The implementation assumes:

• An isotropic \(NVT\) ensemble at equilibrium.

• If per_particle=True (default), the inputs V, U, W, K are per-particle values and if per_particle=False, the inputs are treated as values for the entire system. N is used in the conversion.

• The external temperature (T_ext), \(T\), defaults to \(\langle T_\mathrm{inst}\rangle\) if not provided, where the instantaneous temperature is computed from equipartition, \(T_\mathrm{inst} = \frac{2K}{k_B\, N_d}\) and \(\mathrm{N_d}\) is the number of degrees of freedom (dof).

The isochoric heat capacity (per particle), \(c_v=\frac{1}{N}\left(\frac{\partial E}{\partial T}\right)_V\), and thermal pressure coefficient \(\beta_v = \left(\frac{\partial p}{\partial T}\right)_V\) are computed as

\[\begin{split}c_v &= \frac{\mathrm{Var}(E)}{N\, k_B\, T^2}, \\ \beta_V &= \frac{\mathrm{Cov}(p,E)}{k_B T^2}.\end{split}\]

respectively, where the internal energy is \(E = U + K\) and the instantaneous pressure is \(p = (N k_B T + W) / V\)

Parameters:

• N (int) – Number of particles, \(N\).

• dof (int) – Total number of degrees of freedom of the system, \(N_d\).

• V (float) – Volume of simulation box, \(V\).

• U (Iterable[float]) – Time series of potential energy, \(U\).

• W (Iterable[float]) – Time series of configurational virial, \(W\).

• K (Iterable[float]) – Time series of kinetic energy, \(K\).

• k_B (float, optional) – Boltzmann constant of the unit system (default 1.0), \(k_B\).

• T_ext (float, optional) – External (bath) temperature, \(T\). If None, then \(\langle T_\mathrm{inst}\rangle\) is used.

• per_particle (bool, optional) – If True (default), the V, U, W, K inputs are interpreted as per-particle values. If False, they are treated as extensive.

Returns:

A dictionary with scalar properties and response functions reported as intensive.

• density: \(\rho = N /V\).

• specific_volume: \(v = V / N\).

• potential_energy: Average potential energy per particle, \(\langle U\rangle/N\).

• kinetic_energy: Average kinetic energy per particle, \(\langle K\rangle/N\).

• internal_energy: Average internal energy per particle, \(\langle E\rangle/N\) where \(E = U + K\).

• kinetic_temperature: Instantaneous kinetic temperature, \(\langle T_\textrm{inst}\rangle/N\) where \(T_\textrm{inst} = 2K/k_B\, N_d\).

• external_temperature: External bath temperature (if provided), \(T\).

• pressure: Average instantaneous pressure, \(\langle p\rangle\), where \(p = (N k_B T_\textrm{inst} + W) / V\).

• isochoric_heat_capacity: Isochoric heat capacity per particle, \(c_V = \frac{1}{N}\left(\frac{\partial E}{\partial T}\right)_V = \textrm{Var}(E)/(N k_B T^2)\).

• isochoric_heat_capacity_excess: Excess (configurational) heat capacity per particle, \(c_V^\textrm{ex} = \frac{1}{N}\left(\frac{\partial U}{\partial T}\right)_V = \textrm{Var}(U)/(N k_B T^2)\).

• thermal_pressure_coefficient: Thermal pressure coefficient \(\beta_V = \left(\frac{\partial p}{\partial T}\right)_V = \textrm{Cov}(p,E)/(k_B T^2)\).

• thermal_pressure_coefficient_excess: Excess thermal pressure coefficient \(\beta_V^\textrm{ex} = \frac{1}{V} \left(\frac{\partial W}{\partial T}\right)_V = \textrm{Cov}(W,U)/(V k_B T^2)\).

• thermodynamic_gruneisen_parameter: \(\gamma_G = V\left(\frac{\partial p}{\partial E}\right)_V = \left(\frac{\partial \ln T}{\partial \ln \rho}\right)_S = \dfrac{\beta_V}{\rho \, c_V}\).

• configurational_adiabatic_scaling_exponent: \(\gamma = \left(\dfrac{\partial \ln T}{\partial \ln \rho}\right)_{s_\textrm{ex}} = \textrm{Cov}(W,U)/\textrm{Var}(U)\) where \(s_\textrm{ex}\) is the excess entropy.

• canonical_virial_energy_correlation: Pearson correlation coefficient between \(W\) and \(U\), \(R=\textrm{Cov}(W,U)/\sqrt{\textrm{Var}(W)\textrm{Var}(U)}\).

Return type:

dict

gamdpy.tools.calc_dynamics(trajectory, first_block, qvalues=7.5, overlap_distances=0.3)

Compute dynamical properties from a trajectory in a HDF5 file object.

This function processes blocks of saved configurations to evaluate time‐dependent dynamical observables, the mean square displacement (MSD), the non‐Gaussian parameter (alpha2), the self‐intermediate scattering function (Fs), and the self overlap (Qs)

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. default: 7.5

• overlap_distances (float or array‐like of shape (num_types,), optional) – Criteria for defining self overlap, Qs : particles are counted if they moved less than ‘overlap_distances’ If a single float is provided, it is broadcast to all particle types. Default: 0.3

Returns:

results – Dictionary containing dynamical 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.25, overlap_distances=0.3)
>>> dynamics.keys()
dict_keys(['times', 'msd', 'alpha2', 'qvalues', 'Fs', 'overlap_distances', 'Qs', 'count'])

gamdpy.tools.calculate_molecular_center_of_masses(configuration: Configuration, molecule: str)

Compute molecular center-of-mass positions and masses.

Parameters:

• configuration (Configuration) – Simulation configuration instance containing topology, atomic positions array (conf[‘r’]), atomic masses array (conf[‘m’]), image flags (conf.r_im), and simulation box lengths.

• molecule (str) – Name of the molecule type to process.

Returns:

• positions (ndarray, shape (n_molecules, 3)) – Center-of-mass coordinates of each molecule.

• masses (ndarray, shape (n_molecules,)) – Total mass of each molecule.

Notes

The returned positions are not wrapped according to periodic boundary conditions.

gamdpy.tools.calculate_molecular_velocities(configuration: Configuration, molecule: str)

Compute molecular center-of-mass velocities.

Parameters:

• configuration (Configuration) – Configuration instance containing topology, atomic velocities, and atomic masses.

• molecule (str) – Name of the molecule type to process.

Returns:

velocities – Center-of-mass velocity vectors for each molecule.

Return type:

ndarray, shape (n_molecules, 3)

Notes

Velocities are computed by combining atomic velocities and masses for each molecule.

gamdpy.tools.calculate_molecular_dipoles(configuration: Configuration, atom_charges: ndarray, molecule: str)

Compute molecular dipole moments, centers of mass, and masses.

Parameters:

• configuration (Configuration) – A Configuration instance, containing topology, coordinates.

• atom_charges (array_like of float, shape (n_atoms,)) – Partial charges for each atom in the specified molecule type.

• molecule (str) – Name of the molecule type to process.

Returns:

• dipoles (ndarray, shape (n_molecules, 3)) – Dipole moment vectors for each molecule.

• positions (ndarray, shape (n_molecules, 3)) – Center-of-mass coordinates of each molecule. These are not wrapped according to periodic boundary conditions.

• masses (ndarray, shape (n_molecules,)) – Total mass of each molecule.

Notes

The returned positions are not wrapped according to periodic boundary conditions.

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: float) → callable

Return a function that returns a constant value.

\[f(x) = y_0,\]

Parameters:

value (float) – The value \(y_0\)

Returns:

A function that can be compiled to the device.

Return type:

callable

gamdpy.make_function_ramp(value0: float, x0: float, value1: float, x1: float) → callable

Create a piecewise‐linear “ramp” function.

\[\begin{split}f(x) = \begin{cases} y_0, & x < x_0,\\ y_0 + \dfrac{y_1 - y_0}{x_1 - x_0}\,(x - x_0), & x_0 \le x \le x_1,\\ y_1, & x > x_1, \end{cases}\end{split}\]

Parameters:

• value0 (float) – The value \(y_0\) for \(x < x_0\).

• x0 (float) – The start point \(x_0\) of the linear region.

• value1 (float) – The value \(y_1\) for \(x > x_1\).

• x1 (float) – The end point \(x_1\) of the linear region.

Returns:

A function implementing the above piecewise behavior.

Return type:

callable

gamdpy.make_function_sin(period: float, amplitude: float, offset: float) → callable

Return a function that returns a sin function,

\[f(x) = y_0 + A \sin(2 \pi x / T),\]

with given period (\(T\)), amplitude (\(A\)) and offset (\(y_0\))

Parameters:

• T (float) – The period, \(T\).

• amplitude (float) – The amplitude, \(A\).

• offset (float) – The offset \(y_0\).

Returns:

A function that can be compiled to the device.

Return type:

callable

Miscellaneous

gamdpy.select_gpu()

Select the first available GPU (necessary on some multi-gpu systems).

gamdpy.get_default_sim(num_timeblocks=8)

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() → dict

Return dictionary with flags default compute flags. The boolean value determines whether the corresponding quantity is computed.

gamdpy.conversion_factors(**kwargs)

Return conversion factors from simulation units to common physical units.

The function defines a set of base simulation units by specifying what one unit of length, energy, and mass correspond to in SI (meters, joules, kilograms). From these three base units it derives such as time, pressure, and more. These are returned as a dictionary of multiplicative conversion factors.

Interpretation of the returned dictionary

Each entry cf[key] is a factor that converts a value expressed in simulation units to the unit indicated by key: Example: to give simulation time in ps then use t_ps = t_sim * cf[‘ps’].

param **kwargs:

Exactly one keyword must be provided for each of the three base units: unit_length*, unit_energy*, and unit_mass*. The suffix determines the input unit, e.g.:

• Length: unit_length_in_Angstrom, unit_length_in_nm, …

• Energy: unit_energy_in_kJ_per_mol, unit_energy_in_K, unit_energy_in_eV, …

• Mass: unit_mass_in_u, unit_mass_in_g, …

If no keyword arguments are given, SI base units are assumed: unit_length = 1 m, unit_energy = 1 J, unit_mass = 1 kg.

Special option: get_possible_inputs=True returns a dictionary mapping all accepted input keyword names to their conversion factors to SI.

returns:

Dictionary of conversion factors from simulation units to various units. The dictionary includes (among others):

• length: m, cm, nm, Å, …

• energy: J, kJ/mol, kcal/mol, K, eV, …

• mass: kg, g, u, g/mol, …

• time: s, ps, fs, ns, …

• pressure: Pa, MPa, GPa, bar, atm, …

• …

rtype:

dict

raises KeyError:

If more than one keyword is provided for length, energy, or mass.

raises KeyError:

If any of length, energy, or mass is not specified (and kwargs are not empty).

Examples

>>> from pprint import pprint
>>> from gamdpy.tools import conversion_factors
>>> cf = conversion_factors(unit_length=1.0, unit_energy=1.0, unit_mass=1.0)  # Standard SI units
>>> cf = conversion_factors()  # Standard SI units (same as above)
>>> cf = conversion_factors(unit_length_in_Angstrom=3.4, unit_energy_in_K=120.0, unit_mass_in_u = 39.948) # Argon units
>>> cf = conversion_factors(unit_length_in_Angstrom=1.0, unit_energy_in_kcal_per_mol=1.0, unit_mass_in_u=1.0)  # Molar units (LAMMPS real units)
>>> cf = conversion_factors(unit_length_in_cm=1.0, unit_energy_in_erg=1.0, unit_mass_in_g=1.0)  # centimetre–gram–second (CGS) system
>>> cf = conversion_factors(unit_length_in_nm=1.0, unit_energy_in_kJ_per_mol=1.0, unit_mass_in_u = 1.0)  # Atomistic SI-like units
>>> cf = conversion_factors(unit_length_in_Angstrom=1.0, unit_energy_in_eV=1.0, unit_mass_in_u = 1.0)  # Metallic units

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(num_timeblocks=2)
>>> 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)
restarts/ (Group)
    restart0000/ (Group)
        integrator_state  (Dataset, shape=(1,), dtype=float32)
        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)
    restart0001/ (Group)
        integrator_state  (Dataset, shape=(1,), dtype=float32)
        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)
scalars/ (Group)
    scalars  (Dataset, shape=(2, 64, 3), dtype=float32)
    steps  (Dataset, shape=(65,), dtype=int32)
trajectory/ (Group)
    images  (Dataset, shape=(2, 12, 2048, 3), dtype=int32)
    positions  (Dataset, shape=(2, 12, 2048, 3), dtype=float32)
    ptypes  (Dataset, shape=(2, 12, 2048), dtype=int32)
    steps  (Dataset, shape=(12,), dtype=int32)
    topologies/ (Group)
        block0000/ (Group)
            angles  (Dataset, shape=(0,), dtype=int32)
            bonds  (Dataset, shape=(0,), dtype=int32)
            dihedrals  (Dataset, shape=(0,), dtype=int32)
            molecules/ (Group)
        block0001/ (Group)
            angles  (Dataset, shape=(0,), dtype=int32)
            bonds  (Dataset, shape=(0,), dtype=int32)
            dihedrals  (Dataset, shape=(0,), dtype=int32)
            molecules/ (Group)

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(num_timeblocks=2)
>>> 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 /restarts/:
    - timeblocks_between_restarts: 1
Attributes at /restarts/restart0000/:
    - simbox_data: [12.815602 12.815602 12.815602]
    - simbox_name: Orthorhombic
Attributes at /restarts/restart0000/scalars:
    - scalar_columns: ['U' 'W' 'K' 'm']
Attributes at /restarts/restart0000/topology/molecules/:
    - names: []
Attributes at /restarts/restart0000/vectors:
    - vector_columns: ['r' 'v' 'f']
Attributes at /restarts/restart0001/:
    - simbox_data: [12.815602 12.815602 12.815602]
    - simbox_name: Orthorhombic
Attributes at /restarts/restart0001/scalars:
    - scalar_columns: ['U' 'W' 'K' 'm']
Attributes at /restarts/restart0001/topology/molecules/:
    - names: []
Attributes at /restarts/restart0001/vectors:
    - vector_columns: ['r' 'v' 'f']
Attributes at /scalars/:
    - compression_info: gzip with opts 4
    - scalar_columns: ['U' 'W' 'K']
    - scheduler: Lin
    - scheduler_info: {"steps_between": 16, "npoints": null}
    - steps_between_output: 16
Attributes at /trajectory/:
    - compression_info: gzip with opts 4
    - num_timeblocks: 2
    - scheduler: Log2
    - scheduler_info: {}
    - steps_per_timeblock: 1024
    - trajectory_columns: ['r' 'img']
    - update_ptype: False
    - update_topology: False
Attributes at /trajectory/topologies/block0000/molecules/:
    - names: []
Attributes at /trajectory/topologies/block0001/molecules/:
    - names: []