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)[source]#

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(include_linked_lists=True)[source]#: Autotune the simulation parameters for most efficient calculations on the current machine

run(verbose=True) → None[source]#: Run all blocks of the simulation. See gamdpy.Simulation.run_timeblocks() for an open loop version.

run_timeblocks(num_timeblocks=-1)[source]#

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[source]#

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[source]#: 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'>)[source]#

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[source]#

Scale atom positions and simulation box to a given density.

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

copy_to_device() → None[source]#: Copy all data from host to device memory (CPU to GPU).

copy_to_host() → None[source]#: 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 = False) → Configuration[source]#

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)

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[source]#: Get total potential energy of the configuration

get_volume() → float[source]#: Get volume of simulation box associated with configuration

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

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

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[source]#

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[source]#

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, update_topology: bool = True, verbose: bool = True) → None[source]#

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
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[source]#

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)[source]#

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[source]#: Copy SimulationBox data from host to device.

copy_to_host() → None[source]#: Copy SimulationBox data from device to host.

get_apply_PBC()[source]#: Apply periodic boundary conditions (PBC).

get_dist_moved_exceeds_limit_function()[source]#: For use in neighbor list : Single-particle criterion for whether the neighbor list needs to be rebuilt.

get_dist_sq_dr_function() → callable[source]#: Generates function dist_sq_dr which computes displacement and distance for one neighbor.

get_dist_sq_function() → callable[source]#: Generates function dist_sq_function which computes distance squared for one neighbor.

get_lengths() → ndarray[source]#

Return the box lengths as a numpy array

Returns:: The box lengths in each spatial dimension.
Return type:: numpy.ndarray

get_loop_x_addition()[source]#: 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()[source]#: 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[source]#: Return name of the SimulationBox.

get_volume() → float[source]#

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

Returns:: Volume of the box.
Return type:: float

get_volume_function()[source]#: Returns the function which calculates the volume of the simulation box.

scale(scale_factor: float) → None[source]#: Scale the box lengths by scale_factor

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

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()[source]#: Copy SimulationBox data from host to device.

copy_to_host()[source]#: Copy SimulationBox data from device to host.

get_apply_PBC()[source]#: Apply periodic boundary conditions (PBC).

get_dist_moved_exceeds_limit_function()[source]#: For use in neighbor list : Single-particle criterion for whether the neighbor list needs to be rebuilt.

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

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

get_lengths()[source]#

Return the box lengths as a numpy array

Returns:: The box lengths in each spatial dimension.
Return type:: numpy.ndarray

get_loop_x_addition()[source]#: 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()[source]#: 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()[source]#: Return name of the SimulationBox.

get_volume()[source]#

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

Returns:: Volume of the box.
Return type:: float

get_volume_function()[source]#: Returns the function which calculates the volume of the simulation box.

scale(scale_factor: float) → None[source]#: 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)[source]#

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

Constant temperature and volume#

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

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.

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

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

class gamdpy.Brownian(temperature, tau: float, dt: float, seed=0)[source]#

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)

Constant temperature and pressure#

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

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

class gamdpy.NPT_Langevin(temperature, pressure, alpha: float, alpha_barostat: float, mass_barostat: float, dt: float, volume_velocity=0.0, seed=0)[source]#

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.

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

Other integrators#

class gamdpy.SLLOD(shear_rate, dt)[source]#

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.

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

class gamdpy.GradientDescent(dt: float)[source]#

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

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#

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

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.

Pair potential functions#

gamdpy.LJ_12_6(dist, params)[source]#

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)[source]#

The 12-6 Lennard-Jones potential

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

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

Parameters:

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

gamdpy.harmonic_repulsion(dist, params)[source]#

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)[source]#

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)[source]#

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

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[source]#

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[source]#

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)[source]#

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[source]#

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[source]#

Apply shifted potential cutoff to a pair-potential function

If the input pair potential is $u(r)$, then the shifted potential is $u(r) - u(r_{c})$, where $r_c$ is the cutoff distance. Calls the original potential function twice, avoiding changes to params.

Parameters:: pair_potential (callable) – A function that calculates a pair-potential: u, s, umm = pair_potential(dist, params)
Returns:: pair_potential – A function where shifted 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_force_cutoff(gp.LJ_12_6)
>>> A12, A6, cut = 1.0, 1.0, 2.5
>>> pair_pot = gp.PairPotential(pair_func, params=[A12, A6, cut], max_num_nbs=1000)
>>> interactions = [pair_pot, ]  # List of interactions only containing the pair potential

gamdpy.apply_shifted_force_cutoff(pair_potential)[source]#

Apply shifted force cutoff to a pair-potential function

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

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

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

Fixed interactions#

Classes#

class gamdpy.Bonds(bond_potential, potential_params, indices)[source]#

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

Parameters:

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

See also

gamdpy.harmonic_bond_function: Harmonic bond potential

class gamdpy.Angles(potential, indices, parameters)[source]#

class gamdpy.Tether[source]#

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

Parameters:

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

class gamdpy.Gravity[source]#

Adding a gravitational-like force on particles.

Parameters:

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

class gamdpy.Relaxtemp[source]#

Thermostat using simple kinetic temperature relaxation of each particles.

Parameters:

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

Generators#

gamdpy.make_planar_calculator(configuration, potential_function) → callable[source]#

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[source]#: Returns a dictionary with planar interactions

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

Bond functions#

A bond potential is needed for the Bonds class.

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

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

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)[source]#

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 = 16, compute_flags=None, verbose=False, compression='gzip', compression_opts=4)[source]#: Runtime action for saving scalar data (such as thermodynamic properties) during a timeblock every steps_between_output time steps.

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

class gamdpy.MomentumReset(steps_between_reset: int)[source]#: 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)[source]#: 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)[source]#

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()[source]#

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[source]#

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()[source]#: 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')[source]#

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)[source]#: Generate q-vectors inside a sphere of radius q_max

read(bins) → dict[source]#

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')[source]#: Save average structure factors to a file.

update() → None[source]#: 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')[source]#

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()[source]#

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()[source]#: Update state the current configuration.

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

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)[source]#

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='')[source]#

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[source]#: Returns self.h5

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

load_rumd3(name: str) → File[source]#: 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')[source]#

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[source]#

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[source]#

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) → str[source]#

Convert a configuration to a string formatted as LAMMPS dump file

Parameters:

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

Returns:

string formatted as LAMMPS dump file

Return type:

str

Example

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

Post-analysis tools#

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

Extracts scalar data from simulation output.

Parameters:

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

Returns:

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

Return type:

tuple

Example

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

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

Compute dynamical properties from a saved trajectory HDF5 file.

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

Parameters:

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

Returns:

results – Dictionary containing dynamcal data.

Return type:

dict

Examples

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

Usage within a Python script:

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

gamdpy.tools.calculate_molecular_center_of_masses(configuration: Configuration, molecule: str)[source]#

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)[source]#

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)[source]#

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[source]#

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[source]#

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[source]#

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()[source]#: Select the first available GPU (necessary on some multi-gpu systems).

gamdpy.get_default_sim(num_timeblocks=8)[source]#

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)[source]#

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[source]#: Return dictionary with flags default compute flags. The boolean value determines whether the corresponding quantity is computed.

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

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)[source]#

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)
        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)
        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='/')[source]#

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: []