4. API reference

4.1. Integrator

Integrates equations of motion.

Implements the Velocity Verlet integrator. The reversible reference system propagator algorithm (rRESPA) integrator uses the Velocity Verlet functions inside the MD loop in the main module.

hymd.integrator.integrate_position(positions, velocities, time_step)

Position update step of the Velocity Verlet integration algorithm.

Computes the position update step of the Velocity Verlet algorithm. The positions argument is not changed in place.

Parameters:

positions(N, D) numpy.ndarray

Array of positions of N particles in D dimensions.

velocities(N, D) numpy.ndarray

Array of velocities of N particles in D dimensions.

time_stepfloat

The time step used in the integration.

See also

integrate_velocity

The velocity update step of the Velocity Verlet algorithm.

Notes

The Velocity Verlet algorithm contains two steps, first the velocties are integrated one half step forward in time by applying the forces from the previous step, then the positions are updated a full step using the updated velocities

\[\mathbf{x}_{\text{new}} = \mathbf{x} + \Delta t \mathbf{v}\]

hymd.integrator.integrate_velocity(velocities, accelerations, time_step)

Velocity update step of the Velocity Verlet integration algorithm.

Computes the velocity update step of the Velocity Verlet algorithm. The velocities argument is not changed in place.

Parameters:

velocities(N, D) numpy.ndarray

Array of velocities of N particles in D dimensions.

accelerations(N, D) numpy.ndarray

Array of accelerations of N particles in D dimensions.

time_stepfloat

The time step used in the integration.

See also

integrate_position

The position update step of the Velocity Verlet algorithm.

Notes

The Velocity Verlet algorithm contains two steps, first the velocties are integrated one half step forward in time by applying the forces from the previous step

\[\mathbf{v}_{\text{new}} = \mathbf{v} + \frac{\Delta t}{2m} \mathbf{f}\]

before the positions are moved a full step using the updated half-step velocities.

4.2. Thermostat

Scales or otherwise modifies the particle velocities during simulation to simulate coupling to an external heat bath with temperature T₀.

hymd.thermostat._random_chi_squared(prng: Generator, M: int) → float

Draw the sum of M squared normally distributed values

The value is generated by the Gamma distribution, in lieu of generating M Gaussian distributed numbers and summing their squares.

Parameters:

prngnp.random.Generator

Numpy object that provides a stream of random bits

Mint

Number of standard normally distributed numbers in the sum.

Returns:

float

The sum of M squared normally distributed values centered at zero with unit standard deviation.

Notes

The sum of the squares of k independent standard normal random variables is distributed according to a \(\chi^2\) distribution.

\[\chi^2_1 + \chi_2^2 + \chi_3^2 + \dots + \chi_M^2 \sim \sigma^2\chi^2(k)\]

This is a special case of the \(\Gamma\) distribution, \(\Gamma(k/2, 2)\), and may be generated by

\[\chi^2(k) \sim 2 \Gamma(k/2, 2)\]

References

Knuth, D.E. 1981, Seminumerical Algorithms, 2nd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), pp. 120ff. J. H. Ahrens and U. Dieter, Computing 12 (1974), 223-246.

hymd.thermostat._random_gaussian(prng: Generator) → float

Draw a single random number from the standard normal distribution Generate a number from the Gaussian distribution centered at zero with unit standard deviation, \(N(0, 1)\).

Parameters:

prngnp.random.Generator

Numpy object that provides a stream of random bits

Returns:

float

A random number drawn from \(N(0, 1)\).

hymd.thermostat.csvr_thermostat(velocity: ~numpy.ndarray, names: ~numpy.ndarray, config: ~hymd.input_parser.Config, prng: ~numpy.random._generator.Generator, comm: ~mpi4py.MPI.Intracomm = <mpi4py.MPI.Intracomm object>, random_gaussian: ~typing.Callable[[~numpy.random._generator.Generator], float] = <function _random_gaussian>, random_chi_squared: ~typing.Callable[[~numpy.random._generator.Generator, int, float], float] = <function _random_chi_squared>, remove_center_of_mass_momentum: bool = True) → ndarray

Canonical sampling through velocity rescaling thermostat

Implements the CSVR thermostat. Rescales the system kinetic energy by a stochastically chosen factor to keep the temperature constant. Requires communcation of the kinetic energies calculated locally for each MPI rank. The random numbers sampled through random_gaussian and random_chi_squared are broadcast from the root rank to the other ranks to ensure the scaling is performed with the same stochasticity for all particles in the full system.

The velocities are cleaned of center of mass momentum before the thermostat is applied, and the center of mass momentum is subsequently reapplied. This is performed for each thermostat coupling group, i.e. the center of mass momenta of each group is separately removed and reapplied after thermostatting.

The implementation here is based on the derivation presented in the 2008 Comput. Phys. Commun paper, not in the original 2007 J. Chem. Phys. paper.

Parameters:

velocity(N, D) numpy.ndarray

Array of velocities of N particles in D dimensions.

names(N,) numpy.ndarray

Array of particle names.

configConfig

Configuration dataclass containing simulation metadata and parameters.

prngnp.random.Generator

Numpy object that provides a stream of random bits

commMPI.Intracomm, optional

MPI communicator to use for rank commuication. Defaults to MPI.COMM_WORLD.

random_gaussiancallable, optional

Function for generating standard normally distributed numbers.

random_chi_squaredcallable, optional

Function for generating \(\chi^2\)-distributed numbers

remove_center_of_mass_momentumbool, optional

If True, the center of mass of each coupling group is removed before the thermostat is applied. The center of mass momenta are added back after the kinetic energy rescaling is complete.

See also

_random_gaussian

Used to sample Gaussian-distributed numbers.

_random_chi_squared

Used to sample \(\chi^2\)-distributed numbers.

hymd.input_parser.Config

Configuration dataclass handler.

References

G. Bussi, D. Donadio, and M. Parrinello, J. Chem. Phys. 126, 014101 (2007). G. Bussi and M. Parrinello, Comput. Phys. Commun. 179, 26-29, (2008).

4.3. Hamiltonian

class hymd.hamiltonian.Hamiltonian(config)

Interaction energy functional superclass

__init__(config)

Constructor

Parameters:

configConfig

Configuration object.

See also

hymd.input_parser.Config

Configuration dataclass handler.

_setup()

Superclass setup

Sets up the grid-independent filtering function H, and the SymPy logic for symbolically differentiating interaction energy functionals in Hamiltonian subclasses.

class hymd.hamiltonian.SquaredPhi(config)

Simple squared density interaction energy functional

The interaction energy density takes the form

\[w[\tilde\phi] = \frac{1}{2\kappa\rho_0} \left( \sum_k \tilde\phi_k \right)^2,\]

where \(\kappa\) is the compressibility and \(rho_0\) is the average density of the fully homogenous system. Expressing the species densities in terms of fluctuations from the average,

\[\mathrm{d}\tilde\phi_k = \rho_0 - \tilde\phi_k,\]

it is evident that this interaction energy functional is a slight change of the DefaultNoChi Hamiltonian, as the expanded interaction energy density becomes

\[w[\tilde\phi] = \frac{1}{2\kappa\rho_0} \left( 6\rho_0\left[ \sum_k \mathrm{d}\tilde\phi_k \right] + 9\rho_0^2 + \sum_k \mathrm{d}\tilde\phi_k^2 + \prod_{k\not=l} \mathrm{d}\tilde\phi_k \mathrm{d}\tilde\phi_l \right),\]

identical to DefaultNoChi apart from a constant energy shift \(9\rho_0^2`(which does not impact dynamics) and the :math:\)rho_0`–\(\mathrm{d}\tilde\phi_k\) cross-terms. These cross-terms constitute contributions to the energy linear in \(\mathrm{d}\tilde\phi_k\) absent in DefaultNoChi, which has only quadratic terms present.

See also

hymd.hamiltonian.DefaultNoChi

__init__(config)

Constructor

Parameters:

configConfig

Configuration object.

See also

hymd.input_parser.Config

Configuration dataclass handler.

setup()

Setup the interaction energy potential and the external potential

class hymd.hamiltonian.DefaultNoChi(config)

Incompressibility-only interaction energy functional

The interaction energy density takes the form

\[w[\tilde\phi] = \frac{1}{2\kappa} \left( \sum_k \tilde\phi_k - a \right)^2,\]

where \(\kappa\) is the compressibility and \(a=\rho_0\) for NVT runs where \(\rho_0\) is the average density of the fully homogenous system. In case of NPT runs, \(a\) is a calibrated parameter to obtain the correct average density at the target temperature and pressure. The SquaredPhi Hamiltonian implements a similar functional with an additional linear term component depending on

\[\mathmr{d}\tilde\phi_k = \tilde\phi_k - \rho_0\]

and not \(\mathrm{d}\tilde\phi_k^2\). No explicit inter-species interaction is present apart from the indirect interaction through the incompressibility.

See also

hymd.hamiltonian.DefaultNoChi

hymd.input_parser.Config

Configuration dataclass handler.

__init__(config)

Constructor

Parameters:

configConfig

Configuration object.

See also

hymd.input_parser.Config

Configuration dataclass handler.

setup()

Setup the interaction energy potential and the external potential

class hymd.hamiltonian.DefaultWithChi(config, unique_names, type_to_name_map)

Incompressibility and \(\chi\)-interactions energy functional

The interaction energy density takes the form

\[w[\tilde\phi] = \frac{1}{2\rho_0} \sum_{k,l}\chi_{kl} \tilde\phi_k \tilde\phi_l + \frac{1}{2\kappa} \left( \sum_k \tilde\phi_k - a \right)^2,\]

where \(\kappa\) is the compressibility and \(a=\rho_0\) for NVT runs where \(\rho_0\) is the average density of the fully homogenous system. In case of NPT runs, \(a\) is a calibrated parameter to obtain the correct average density at the target temperature and pressure. \(\chi_{ij}\) is the Flory-Huggins-like inter-species mixing energy.

__init__(config, unique_names, type_to_name_map)

Constructor

Parameters:

configConfig

Configuration object.

unique_namesnumpy.ndarray

Sorted array of all distinct names of different species present in the simulation. Result of numpy.unique(all_names), where all_names is the gathered union of all individual MPI ranks’ names arrays.

type_to_name_mapdict[int, str]

Dictionary of the mapping from type indices (integers) to type names.

See also

hymd.input_parser.Config

Configuration dataclass handler.

setup(unique_names, type_to_name_map)

Setup the interaction energy potential and the external potential

Parameters:

unique_namesnumpy.ndarray

Sorted array of all distinct names of different species present in the simulation. Result of numpy.unique(all_names), where all_names is the gathered union of all individual MPI ranks’ names arrays.

type_to_name_mapdict[int, str]

Dictionary of the mapping from type indices (integers) to type names.

4.4. Force

Calculates intramolecular forces between bonded particles in molecules

class hymd.force.Angle(atom_1: str, atom_2: str, equilibrium: float, strength: float, atom_3: str)

Dataclass representing a single three-particle bond type

A bond type is a bond strength and equilibrium angle associated with any three-particle bond between particles of specific types A, B, and C (where A, B, and C may be different or the same). Harmonic angular three-particle bonds in HyMD take the form

\[V_3(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = \frac{1}{2}k \left( \cos^{-1} \left[ \frac{ (\mathbf{r}_1-\mathbf{r}_2) \cdot (\mathbf{r}_3-\mathbf{r}_2) } { \vert\mathbf{r}_1-\mathbf{r}_2\vert \vert\mathbf{r}_3-\mathbf{r}_2\vert } \right] - \theta_0 \right)^2\]

See also

Bond

Two-particle bond type dataclass

Attributes:

atom_1str

Type name of particle 1.

atom_2str

Type name of particle 2.

atom_3str

Type name of particle 3.

equilibriumfloat

Equilibrium angle at which the energy associated with the three-particle angular bond vanishes and the resulting force is zero.

strengthfloat

Harmonic bond strength coefficient.

atom_3: str

class hymd.force.Bond(atom_1: str, atom_2: str, equilibrium: float, strength: float)

Dataclass representing a single two-particle bond type

A bond type is a bond strength and equilibrium distance associated with any bond between particles of specific types A and B (where A and B may be the same or different). Harmonic two-particle bonds in HyMD take the form

\[V_2(\mathbf{r}_1, \mathbf{r}_2) = \frac{1}{2}k \left( \vert \mathbf{r}_1 - \mathbf{r}_2 \vert - r_0 \right)^2,\]

where \(k\) is the bond strength (spring constant) and \(r_0\) is the equilibrium bond length (at which the energy is zero).

Attributes:

atom_1str

Type name of particle 1.

atom_2str

Type name of particle 2.

equilibriumfloat

Equilibrium distance at which the energy associated with the bond vanishes and the resulting force is zero.

strengthfloat

Harmonic bond strength coefficient (spring constant).

atom_1: str

atom_2: str

equilibrium: float

strength: float

class hymd.force.Chi(atom_1: str, atom_2: str, interaction_energy: float)

Dataclass representing a single \(\chi\) mixing interaction type

An interaction mixing energy type is a mixing energy associated with density overlap between species of types A and B (specified as inputs atom_1 and atom_2). A positive mixing energy promotes phase separation, a negative mixing energy promotes mixing. The interaction energy density (provided the DefaultWithChi Hamiltonian is used) takes the form

\[w[\tilde\phi] = \frac{1}{2\kappa} \sum_{k,l}\chi_{k,l} \tilde\phi_k\tilde\phi_l,\]

where \(\chi_{k,l}\) denotes the mixing energy between species \(k\) and \(l\), with \(\kappa\) being the incompressibility. The value of the interaction mixing energy parameter may be extracted from Flory-Huggins theory.

See also

hymd.hamiltonian.DefaultWithChi

Interaction energy functional using \(\chi\)-interactions.

Attributes:

atom_1str

Type name of particle 1.

atom_2str

Type name of particle 2.

interaction_energyfloat

Interaction mixing energy.

atom_1: str

atom_2: str

interaction_energy: float

class hymd.force.Dielectric_type(atom_1: str, dielectric_value: float)

atom_1: str

dielectric_value: float

class hymd.force.Dihedral(atom_1: str, atom_2: str, atom_3: str, atom_4: str, coeffs: ndarray, dih_type: int)

Dataclass representing a single four-particle bond type

A bond type is a bond strength and equilibrium angle associated with any four-particle torsional bond between particles of specific types A, B, C, and D (where A, B, C, and D may be different or the same). Dihedral four-particle bonds in HyMD take different forms depending on the dih_type parameter.

In the following, let \(\phi\) denote the angle between the planes spanned by the relative positions of atoms A-B-C and B-C-D. If dih_type = 0, then

\[V_4(\phi) = \sum_n c_n \left( 1 + \cos\left[ n\phi - d_n \right] \right),\]

where \(c_n\) are energy coefficients and \(d_n\) are propensity phases. By default, the cosine sum is truncated at five terms. If coeffs provides only a single float, this is used as the coiling propensity parameter \(\lambda\). In this case, \(c_n\) and \(d_n\) are automatically set to values which promote alpha helical (\(\lambda=-1\)), beta sheet (\(\lambda=1\)), or a mixed (\(1>\lambda>-1\)) structure (using values provided by Bore et al. (2018)). In this case, the full potential takes the form

\[V_4(\phi;\lambda) = \frac{1}{2}(1-\lambda) V_{\text{prop},\alpha}(\phi) + \frac{1}{2}(1+\lambda) V_{\text{prop},\beta}(\phi) + (1-\vert\lambda\vert) V_{\text{prop}, \text{coil}}(\phi),\]

with each of the \(V_{\mathrm{prop}, X}\) being fixed cosine series potentials with pre-set \(c_n\) -s and \(d_n\) -s.

If dih_type = 1, then a combined bending torsional (CBT) potential is employed,

\[V_4(\phi,\gamma;\lambda) = V_4(\phi;\lambda) + \frac{1}{2}K(\phi)(\gamma - \gamma_0)^2,\]

where \(V_4(\phi;\lambda)\) specifies the potential as given by the coiling propensity parameter above, \(K(\phi)\) is a separate cosine series potential acting as the angular three-particle bond strength, and \(\gamma_0\) is the three-particle bond equilibrium angle. In this case, the coeffs parameter must specify both a \(\lambda\) value, in additon to the energy and phases dictating the \(K(\phi)\) potential.

References

Bore et al. J. Chem. Theory Comput., 14(2): 1120–1130, 2018.

Attributes:

atom_1str

Type name of particle 1.

atom_2str

Type name of particle 2.

atom_3str

Type name of particle 3.

atom_4str

Type name of particle 4.

coeffslist[list[float]] or list[float] or numpy.ndarray

Dihedral coefficients defining the series expansion of the dihedral energy.

dih_typeint

Specifies the type of dihedral used; If 0, then coeffs must contain either two lists of five energy coefficients \(c_n\) and five propensity phases \(d_n\) or a single floating point number defining the \(\lambda\) coiling propensity parameter. If 1, the combined bending-torsional potential is used, and coeffs must specify \(\lambda\) and two lists containing energy coefficients and propensity phases for the \(V_\text{prop}\) propensity potential, defined in terms of cosine series.

atom_1: str

atom_2: str

atom_3: str

atom_4: str

coeffs: ndarray

dih_type: int

hymd.force.compute_angle_forces__plain(f_angles, r, bonds_3, box_size)

Computes forces resulting from angular interactions

Deprecated since version 1.0.0: compute_angle_forces__plain was replaced by compiled Fortran code prior to 1.0.0 release.

hymd.force.compute_bond_forces__plain(f_bonds, r, bonds_2, box_size)

Computes forces resulting from bonded interactions

Deprecated since version 1.0.0: compute_bond_forces__plain was replaced by compiled Fortran code prior to 1.0.0 release.

hymd.force.compute_dihedral_forces__plain(f_dihedrals, r, bonds_4, box_size)

Computes forces resulting from dihedral interactions

Deprecated since version 1.0.0: compute_dihedral_forces__plain was replaced by compiled Fortran code prior to 1.0.0 release.

hymd.force.dipole_forces_redistribution(f_on_bead, f_dipoles, trans_matrices, a, b, c, d, type_array, last_bb)

Redistribute electrostatic forces calculated from topologically reconstructed ghost dipole point charges to the backcone atoms of the protein.

hymd.force.find_all_paths(G, u, n)

Helper function that recursively finds all paths of given lenght ‘n + 1’ inside a network ‘G’. Adapted from https://stackoverflow.com/a/28103735.

hymd.force.prepare_bonds(molecules, names, bonds, indices, config, topol=None)

Rearrange the bond information for usage in compiled Fortran kernels

Restructures the lists resulting from the execution of prepare_bonds_old into numpy arrays suitable for calls to optimized Fortran code calculating bonded forces and energiesself.

Parameters:

molecules(N,) numpy.ndarray

Array of integer molecule affiliation for each of N particles. Global (across all MPI ranks) or local (local indices on this MPI rank only) may be used, both, without affecting the result.

names(N,) numpy.ndarray

Array of type names for each of N particles.

bonds(N,M) numpy.ndarray

Array of M bonds originating from each of N particles.

indices(N,) numpy.ndarray

Array of integer indices for each of N particles. Global (across all MPI ranks) or local (local indices on this MPI rank only) may be used, both, without affecting the result.

configConfig

Configuration object.

Returns:

bonds_2_atom_1(B,) numpy.ndarray

Local index of particle 1 for each of B constructed two-particle bonds.

bonds_2_atom_2(B,) numpy.ndarray

Local index of particle 2 for each of B constructed two-particle bonds.

bonds_2_equilibrium(B,) numpy.ndarray

Equilibrium distance for each of B constructed two-particle bonds.

bonds_2_strength(B,) numpy.ndarray

Bond strength for each of B constructed two-particle bonds.

bonds_3_atom_1(A,) numpy.ndarray

Local index of particle 1 for each of A constructed three-particle bonds.

bonds_3_atom_2(A,) numpy.ndarray

Local index of particle 2 for each of A constructed three-particle bonds.

bonds_3_atom_3(A,) numpy.ndarray

Local index of particle 3 for each of A constructed three-particle bonds.

bonds_3_equilibrium(A,) numpy.ndarray

Equilibrium angle for each of A constructed three-particle bonds.

bonds_3_strength(A,) numpy.ndarray

Bond strength for each of A constructed three-particle bonds.

bonds_4_atom_1(D,) numpy.ndarray

Local index of particle 1 for each of D constructed four-particle bonds.

bonds_4_atom_2(D,) numpy.ndarray

Local index of particle 2 for each of D constructed four-particle bonds.

bonds_4_atom_3(D,) numpy.ndarray

Local index of particle 3 for each of D constructed four-particle bonds.

bonds_4_atom_4(D,) numpy.ndarray

Local index of particle 4 for each of D constructed four-particle bonds.

bonds_4_coeff(D,) numpy.ndarray

Cosine series coefficients for each of D constructed four-particle bonds.

bonds_4_type(D,) numpy.ndarray

Dihedral type for each of D constructed four-particle bonds.

bonds_4_last(D,) numpy.ndarray

Flags indicating if dih_type is 1 for each of D constructed four-particle bonds.

See also

prepare_bonds_old

Used internally to reconstruct the bonded interactions types from the connectivity information in the structure/topology input file and the bonded types specified in the configuration file.

hymd.force.prepare_bonds_old(molecules, names, bonds, indices, config)

Find bonded interactions from connectivity and bond types information

Deprecated since version 1.0.0: prepare_bonds_old was replaced by prepare_bonds for use with compiled Fortran kernels prior to 1.0.0 release.

Prepares the necessary equilibrium and bond strength information needed by the intramolecular interaction functions. This is performed locally on each MPI rank, as the domain decomposition ensures that for all molecules all consituent particles are always contained on the same MPI rankself.

This function traverses the bond connectivity information provided in the structure/topology input file and indentifies any two-, three-, or four-particle potential bonds. For each connected chain of two, three, or four particles, a matching to bond types is attempted. If the corresponding names match, a bond object is initialized.

In order to investigate the connectivity, a graph is created using networkx functionality.

Parameters:

molecules(N,) numpy.ndarray

Array of integer molecule affiliation for each of N particles. Global (across all MPI ranks) or local (local indices on this MPI rank only) may be used, both, without affecting the result.

names(N,) numpy.ndarray

Array of type names for each of N particles.

bonds(N,M) numpy.ndarray

Array of M bonds originating from each of N particles.

indices(N,) numpy.ndarray

Array of integer indices for each of N particles. Global (across all MPI ranks) or local (local indices on this MPI rank only) may be used, both, without affecting the result.

configConfig

Configuration object.

Returns:

bonds_2list

List of lists containing local particle indices, equilibrium distance, and bond strength coefficient for each reconstructed two-particle bond.

bonds_3

List of lists containing local particle indices, equilibrium angle, and bond strength coefficient for each reconstructed three-particle bond.

bonds_4

List of lists containing local particle indices, dihedral type index, and dihedral coefficients for each reconstructed four-particle torsional bond.

bb_indexlist

List indicating the dihedral type of each four-particle bond in bonds_4.

See also

Bond

Two-particle bond type dataclass.

Angle

Three-particle angular bond type dataclass.

Dihedral

Four-particle torsional bond type dataclass.

hymd.input_parser.Config

Configuration dataclass handler.

hymd.force.prepare_index_based_bonds(molecules, topol)

hymd.force.propensity_potential_coeffs(x: float)

4.5. Logger

class hymd.logger.Logger

Log output handler class

Notes

This wraps the default python library logging, see docs.python.org/3/library/logging.html.

Attributes:

levelint

Determines the verbosity level of the log, corresponding to logging.level. Numerical values 50 (logging.CRITICAL), 40 (logging.ERROR), 30 (logging.WARNING), 20 (logging.INFO), 10 (logging.DEBUG), and 0 (logging.UNSET) are supported values. Any log event message less severe than the specified level is ignored. All other event messages are emitted.

log_filestr

Path to output log file.

formatstr

Prepended dump string for each log event. Specifies the log event level, the module emitting the event, the code line, the enclosing function name, and the MPI rank writing the message.

date_formatstr

Prepends the date before all log event messages.

formatterlogging.Formatter

Formatter handling the prepending of the information in format and date_format to each log event message. Used by default for all loggers.

rank0logging.Logger

Default logger object for the root MPI rank.

all_rankslogging.Logger

Default logger object for messages being emitted from all MPI ranks simultaneously.

classmethod setup(default_level=20, log_file=None, verbose=False, comm=<mpi4py.MPI.Intracomm object>)

Sets up the logger object.

If a log_file path is provided, log event messages are output to it. Otherwise, the logging messages are emitted to stdout.

Parameters:

default_levelint, optional

Default verbosity level of the logger. Unless specified, it is 10 (logging.INFO).

log_filestr, optional

Path to output log file. If None or not priovided, no log file is used and all logging is done to stdout.

verbosebool, optional

Increases the logging level to 30 (logging.WARNING) if True, otherwise leaves the logging level unchanged.

class hymd.logger.MPIFilterRoot(comm=<mpi4py.MPI.Intracomm object>)

Log output Filter wrapper class for the root MPI rank log

filter(record)

Log event message filter

Parameters:

recordlogging.LogRecord

LogRecord object corresponding to the log event.

class hymd.logger.MPIFilterAll(comm=<mpi4py.MPI.Intracomm object>)

Log output Filter wrapper class for the all-MPI-ranks log

filter(record)

Log event message filter

Parameters:

recordlogging.LogRecord

LogRecord object corresponding to the log event.

4.6. Input parser

Parses and handles the configuration information provided for the simulation

class hymd.input_parser.Config(n_steps: int, time_step: float, mesh_size: ~typing.List[int] | ~numpy.ndarray | int, sigma: float, kappa: float, dtype: ~numpy.dtype | None = None, box_size: ~typing.List[float] | ~numpy.ndarray | None = None, n_print: int | None = None, tau: float | None = None, start_temperature: float | bool | None = None, target_temperature: float | bool | None = None, mass: float | None = None, hamiltonian: str | None = None, domain_decomposition: int | bool | None = None, integrator: str | None = None, respa_inner: int = 1, file_name: str = '<config file path unknown>', name: str | None = None, tags: ~typing.List[str] = <factory>, bonds: ~typing.List[~hymd.force.Bond] = <factory>, angle_bonds: ~typing.List[~hymd.force.Angle] = <factory>, dihedrals: ~typing.List[~hymd.force.Dihedral] = <factory>, chi: ~typing.List[~hymd.force.Chi] = <factory>, n_particles: int | None = None, max_molecule_size: int | None = None, n_flush: int | None = None, thermostat_work: float = 0.0, thermostat_coupling_groups: ~typing.List[~typing.List[str]] = <factory>, initial_energy: float | None = None, cancel_com_momentum: int | bool = False, coulombtype: str | None = None, convergence_type: str | None = None, pol_mixing: float | None = None, dielectric_const: float | None = None, conv_crit: float | None = None, dielectric_type: ~typing.List[~hymd.force.Dielectric_type] = <factory>, self_energy: float | None = None, type_charges: ~typing.List[float] | ~numpy.ndarray | None = None, rho0: float | None = None, a: float | None = None, pressure: bool = False, barostat: str | None = None, barostat_type: str | None = None, tau_p: float | None = None, target_pressure: ~typing.List[~hymd.barostat.Target_pressure] = <factory>, n_b: int | None = None, m: ~typing.List[float] = <factory>)

Configuration object

Handles and verifies the simulation configuration specified in the configuration file.

See also

hymd.input_parser.Bond

Two-particle bond type dataclass.

hymd.input_parser.Angle

Three-particle bond type dataclass.

hymd.input_parser.Dihedral

Four-particle bond type dataclass.

Attributes:

gas_constantfloat

Constant value of the gas constant, R (equivalently the Boltzmann constant) in the units used internally in HyMD.

coulomb_constantfloat

Constant value of the Coulomb constant which converts electric field values to forces and electric potential values to energies in the units used internally in HyMD.

n_steps: int

Number of time steps in the simulation.

time_step: float

Outer time step used in the simulation. If the rRESPA intergrator is used, the inner time step (the time step used in the integration of intramolecular bending, stretching, and torsional forces) is time_step / respa_inner.

box_sizelist[float] or (D,) numpy.ndarray

Simulation box size of simulation in D dimensions in units of nanometers.

mesh_sizelist[int] or int or numpy.ndarray

Number of grid points used for the discrete density grid.

sigmafloat

Filter width representing the effective coarse-graining level of the particles in the simulation.

kappafloat

Compressibility parameter used in the relaxed incompressibility term in the interaction energy functional.

n_printint, optional

Frequency of trajectory/energy output to the H5MD trajectory/energy output file (in units of number of time steps).

taufloat, optional

The time scale of the CSVR thermostat coupling.

start_temperaturefloat, optional

Generate starting temperature by assigning all particle velocities randomly according to the Maxwell-Boltzmann distribution at start_temperature Kelvin prior to starting the simulation.

target_temperaturefloat, optional

Couple the system to a heat bath at target_temperature Kelvin.

massfloat, optional

Mass of the particles in the simulation.

hamiltonianstr, optional

Specifies the interaction energy functional \(W[\tilde\phi]\) for use with the particle-field interactions. Options: SquaredPhi, DefaultNohChi, or DefaultWithChi.

domain_decompositionint, optional

Specifies the interval (in time steps) of domain decomposition exchange, involving all MPI ranks sending and receiving particles according to the particles’ positions in the integration box and the MPI ranks’ assigned domain.

integratorstr, optional

Specifies the time integrator used in the simulation. Options: velocity-verlet or respa.

respa_innerint, optional

The number of inner time steps in the rRESPA integrator. This denotes the number of intramolecular force calculations (stretching, bending, torsional) are performed between each impulse applied from the field forces.

file_namestr, optional

File path of the parsed configuration file.

namestr, optional

Name for the simulation.

tagslist[str], optional

Tags for the simulation.

bondslist[Bond], optional

Specifies harmonic stretching potentials between particles in the same molecule.

angle_bondslist[Angle], optional

Specifies harmonic angular bending potentials between particles in the same molecule.

dihedralslist[Dihedral], optional

Specifies four-particle torsional potentials by cosine series.

chilist[Chi], optional

Specifies \(\chi\)-interaction parameters between particle species.

n_particlesint, optional

Specifies the total number of particles in the input. Optional keyword for validation, ensuring the input HDF5 topology has the correct number of particles and molecules.

max_molecule_sizeint, optional

Maximum size of any single molecule in the system. Used to speed up distribution of particles onto MPI ranks in a parallel fashion.

n_flushint, optional

Frequency of HDF5 write buffer flush, forcing trajectory/energy to be written to disk (in units of number of n_print).

thermostat_workfloat

Work performed by the thermostat on the system.

thermostat_coupling_groupslist[str], optional

Specifies individual groups coupling independently to the CSVR thermostat. E.g. in a system containing "A", "B", and "C" type particles, thermostat_coupling_groups = [["A", "B"], ["C"],] would thermalise types "A" and "B" together and couple "C" type particles to a different thermostat (all individual thermostats are at the same temperature, i.e. target_temperature Kelvin).

initial_energyfloat

Value of the total energy prior to the start of the simulation.

cancel_com_momentumint, optional

If True, the total linear momentum of the center of mass is removed before starting the simulation. If an integer is specifed, the total linear momentum of the center of mass is removed every remove_com_momentum time steps. If False, the linear momentum is never removed.

coulombtypestr, optional

Specifies the type of electrostatic Coulomb interactions in the system. The strength of the electrostatic forces is modulated by the relative dielectric constant of the simulation medium, specified with the dielectric_const keyword. Charges for individual particles are specified in the structure/topology HDF5 input file, not in the configuration file. If no charges (or peptide backbone dipoles) are present, the electrostatic forces will not be calculated even if this keyword is set to "PIC_Spectral".

dielectric_constfloat, optional

Specifies the relative dielectric constant of the simulation medium which regulates the strength of the electrostatic interactions. When using helical propensity dihedrals, this keyword must be specified—even if electrostatics are not included with the coulombtype keyword.

dielectric_type: list[float], optional

Specifies the relative dielectric constant of the simulation medium which regulates the strength of the electrostatic interactions. The list assigns relative dielectric values to each bead type, and an anisotropic dielectric function is obtained from a weighted average.

a: float = None

angle_bonds: List[Angle]

barostat: str = None

barostat_type: str = None

bonds: List[Bond]

box_size: List[float] | ndarray = None

cancel_com_momentum: int | bool = False

chi: List[Chi]

conv_crit: float = None

convergence_type: str = None

coulomb_constant: ClassVar[float] = 138.935458

coulombtype: str = None

dielectric_const: float = None

dielectric_type: List[Dielectric_type]

dihedrals: List[Dihedral]

domain_decomposition: int | bool = None

dtype: dtype = None

file_name: str = '<config file path unknown>'

gas_constant: ClassVar[float] = 0.0083144621

hamiltonian: str = None

initial_energy: float = None

integrator: str = None

kappa: float

m: List[float]

mass: float = None

max_molecule_size: int = None

mesh_size: List[int] | ndarray | int

n_b: int = None

n_flush: int = None

n_particles: int = None

n_print: int = None

n_steps: int

name: str = None

pol_mixing: float = None

pressure: bool = False

respa_inner: int = 1

rho0: float = None

self_energy: float = None

sigma: float

start_temperature: float | bool = None

tags: List[str]

target_pressure: List[Target_pressure]

target_temperature: float | bool = None

tau: float = None

tau_p: float = None

thermostat_coupling_groups: List[List[str]]

thermostat_work: float = 0.0

time_step: float

type_charges: List[float] | ndarray = None

hymd.input_parser.check_NPT_conditions(config, comm=<mpi4py.MPI.Intracomm object>)

Check validity of barostat_type, barostat, a, rho0, target_pressure, tau_p

hymd.input_parser.check_angles(config, names, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_bonds(config, names, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_box_size(config, input_box, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_cancel_com_momentum(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_charges(config, charges, comm=<mpi4py.MPI.Intracomm object>)

Check if charges across ranks sum to zero.

Parameters:

charges(N,) numpy.ndarray

Array of floats with charges for N particles.

commmpi4py.Comm, optional

MPI communicator, defaults to mpi4py.COMM_WORLD.

hymd.input_parser.check_charges_types_list(config, types, charges, comm=<mpi4py.MPI.Intracomm object>)

Creates a list of charge values of length types. Charges are sorted according to type ID. Used in field.py. # TODO: this is messy, we should fix it

hymd.input_parser.check_chi(config, names, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_config(config, indices, names, types, charges, input_box, comm=<mpi4py.MPI.Intracomm object>)

Performs various checks on the specfied config to ensure consistency

Parameters:

configConfig

Configuration object.

indices(N,) numpy.ndarray

Array of integer indices for N particles.

names(N,) numpy.ndarray

Array of string names for N particles.

types(N,) numpy.ndarray

Array of integer type indices for N particles.

charges(N,) numpy.ndarray

Array of floats charges for N particles.

commmpi4py.Comm, optional

MPI communicator, defaults to mpi4py.COMM_WORLD.

Returns:

configConfig

Validated configuration object.

hymd.input_parser.check_dielectric(config, comm=<mpi4py.MPI.Intracomm object>)

Error handling for electrostatics. Unit testing of toml/tomli input.

hymd.input_parser.check_dihedrals(config, names, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_domain_decomposition(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_hamiltonian(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_integrator(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_m(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_mass(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_max_molecule_size(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_n_b(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_n_flush(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_n_particles(config, indices, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_n_print(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_name(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_start_and_target_temperature(config, comm=<mpi4py.MPI.Intracomm object>)

Validate provided starting and target thermostat temperatures

Assesses the provided temperature target and ensures it is a non-negative floating point number or False. Ensures the starting temperature is a non-negative floating point number or False.

If the value for either is None, the returned configuration object has the values defaulted to False in each case.

Parameters:

configConfig

Configuration object.

Returns:

validated_configConfig

Configuration object with validated target_temperature and start_temperature.

hymd.input_parser.check_tau(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.check_thermostat_coupling_groups(config, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.parse_config_toml(toml_content, file_path=None, comm=<mpi4py.MPI.Intracomm object>)

hymd.input_parser.read_config_toml(file_path)

hymd.input_parser.sort_dielectric_by_type_id(config, charges, types)

Creates a list of length N CG-particles, sorted after the charges from the input HDF5 file. Used in file_io.py

4.7. Field

Forces and energies from the discretized particle-field grid interactions

hymd.field.compute_field_and_kinetic_energy(phi, phi_q, psi, velocity, hamiltonian, positions, types, v_ext, config, layouts, comm=<mpi4py.MPI.Intracomm object>)

Compute the particle-field and kinetic energy contributions

Calculates the kinetic energy through

\[E_k = \sum_j \frac{1}{2}m_j\mathbf{v}_j\cdot\mathbf{v}_j,\]

and the particle-field energy by

\[E_\text{field} = \int\mathrm{d}\mathbf{r}\, w[\tilde\phi],\]

where \(w\) is the interaction energy functional density.

Parameters:

philist[pmesh.pm.RealField]

Pmesh RealField objects containing discretized particle number density values on the computational grid; one for each particle type. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

velocity(N,D) numpy.ndarray

Array of velocities for N particles in D dimensions. Local for each MPI rank.

hamiltonianHamiltonian

Particle-field interaction energy handler object. Defines the grid-independent filtering function, \(H\).

positions(N,D) numpy.ndarray

Array of positions for N particles in D dimensions. Local for each MPI rank.

types(N,) numpy.ndarray

Array of type indices for each of N particles. Local for each MPI rank.

v_extlist[pmesh.pm.RealField]

Pmesh RealField objects containing discretized particle-field external potential values on the computational grid; one for each particle type. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

configConfig

Configuration object.

layoutslist[pmesh.domain.Layout]

Pmesh communication layout objects for domain decompositions of each particle type. Used as blueprint by pmesh.pm.readout for exchange of particle information across MPI ranks as necessary.

commmpi4py.Comm

MPI communicator to use for rank commuication.

See also

update_field

Computes the up-to-date external potential for use in calculating the particle-field energy.

hymd.field.compute_field_energy_q_GPE(config, phi_eps, field_q_energy, dot_elec, comm=<mpi4py.MPI.Intracomm object>)

Compute the electrostatic energy after electrosatic forces is calculated.

From the definition of the elecrostatic potential \(\Psi\), the energy is

W = frac{1}{2}intmathrm{d}mathbf{r},

epsilon(mathbf{r})} left(mathbf{E}cdot mathbf{E}right),

where \(\epsilon(\mathbf{r})}\) is the anisotropic, spatially dependent, relative dielectric of the simulation medium.

Parameters:

confighymd.input_parser.Config

Configuration object.

phi_epspmesh.pm.RealField

Pmesh RealField object for storing calculated discretized relative dielectric values on the computational grid. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

field_q_energyfloat

Total elecrostatic energy.

dot_elecpmesh.pm.RealField

Pmesh RealField object for storing \(|\mathbf{E(r)}|^{2}\) on the computational grid. Local for each MPI rank – the full computational grid is represented by the collective fields of all MPI ranks.

commmpi4py.Comm

MPI communicator to use for rank commuication.

See also

update_field_force_q_GPE

Compute the electrosatic force from an anisotropic dielectric general Poisson equation.

hymd.field.compute_field_force(layouts, r, force_mesh, force, types, n_types)

Interpolate particle-field forces from the grid onto particle positions

Backmaps the forces calculated on the grid to particle positions using the window function \(P\) (by default cloud-in-cell [CIC]). In the following, let \(\mathbf{F}_j\) denote the force acting on the particle with index \(j\) and position \(\mathbf{r}_j\). The interpolated force is

\[\mathbf{F}_j = -\sum_k\nabla V_{j_k} P(\mathbf{r}_{j_k}-\mathbf{r}_j)h^3,\]

where \(V_{j_k}\) is the discretized external potential at grid vertex \(j_k\), \(\mathbf{r}_{j_k}\) is the position of the grid vertex with index \(j_k\), and \(h^3\) is the volume of each grid voxel. The sum is taken over all closest neighbour vertices, \(j_k\).

Parameters:

layoutslist[pmesh.domain.Layout]

Pmesh communication layout objects for domain decompositions of each particle type. Used as blueprint by pmesh.pm.readout for exchange of particle information across MPI ranks as necessary.

r(N,D) numpy.ndarray

Array of positions for N particles in D dimensions. Local for each MPI rank.

force_meshlist[pmesh.pm.RealField]

Pmesh RealField objects containing discretized particle-field force density values on the computational grid; D fields in D dimensions for each particle type. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

force(N,D) numpy.ndarray

Array of forces for N particles in D dimensions. Local for each MPI rank.

types(N,) numpy.ndarray

Array of type indices for each of N particles. Local for each MPI rank.

n_typesint

Number of different unique types present in the simulation system. n_types is global, i.e. the same for all MPI ranks even if some ranks contain zero particles of certain types.

hymd.field.compute_self_energy_q(config, charges, comm=<mpi4py.MPI.Intracomm object>)

Compute the self energy for the interaction due to the Ewald scheme used to compute the electrostatics. The energy is stored in config.

The self interaction energy is given by:

\[U_{self} = \sqrt{\frac{1}{2\pi\sigma^2}} \sum_{i=1}^{N} q_i^2\]

where \(q_i\) are the charges and \(\sigma\) is the half-width of the Gaussian filter.

Parameters:

configConfig

Configuration object.

charges(N,) numpy.ndarray

Array of particle charge values for N particles. Local for each MPI rank.

commmpi4py.Comm

MPI communicator to use for rank communication.

Returns:

field_q_self_energyfloat

Electrostatic self energy.

hymd.field.domain_decomposition(positions, pm, *args, molecules=None, bonds=None, topol=False, verbose=0, comm=<mpi4py.MPI.Intracomm object>)

Performs domain decomposition

Rearranges the MPI rank memory layout into one that better mimics the PFFT pencil grid 2D decomposition. As molecules are always required to be fully contained on a single MPI rank, a perfect decomposition is not always possible. The best decomposition which leaves the center of mass of all molecules in their respective correct domains is used, with possible overlap into neighbouring domains if the spatial extent of any molecule crosses domain boundaries.

Parameters:

positions(N,D) numpy.ndarray

Array of positions for N particles in D dimensions. Local for each MPI rank.

pmpmesh.pm.ParticleMesh

Pmesh ParticleMesh object interfacing to the CIC window function and the PFFT discrete Fourier transform library.

*args

Variable length argument list containing arrays to include in the domain decomposition.

molecules(N,) numpy.ndarray, optional

Array of integer molecule affiliation for each of N particles. Global (across all MPI ranks) or local (local indices on this MPI rank only) may be used, both, without affecting the result.

bonds(N,M) numpy.ndarray, optional

Array of M bonds originating from each of N particles.

verboseint, optional

Specify the logging event verbosity of this function.

commmpi4py.Comm

MPI communicator to use for rank commuication.

Returns:

Domain decomposed positions and any array specified in

*args, in addition to bonds and molecules if these

arrays were provided as input.

hymd.field.initialize_pm(pmesh, config, comm=<mpi4py.MPI.Intracomm object>)

Creates the necessary pmesh objects for pfft operations.

Parameters:

pmeshmodule ‘pmesh.pm’

configConfig

Configuration dataclass containing simulation metadata and parameters.

commMPI.Intracomm, optional

MPI communicator to use for rank commuication. Defaults to MPI.COMM_WORLD.

Returns:

pmobject ‘pmesh.pm.ParticleMesh’

field_listlist[pmesh.pm.RealField], (multiple)

Essential list of pmesh objects required for MD

list_coulomblist[pmesh.pm.RealField], (multiple)

Additional list of pmesh objects required for electrostatics.

hymd.field.update_field(phi, phi_laplacian, phi_transfer, layouts, force_mesh, hamiltonian, pm, positions, types, config, v_ext, phi_fourier, v_ext_fourier, m, compute_potential=False)

Calculate the particle-field potential and force density

If compute_potential is True, the energy may subsequently be computed by calling compute_field_and_kinetic_energy.

Computes the particle-field external potential \(V_\text{ext}\) from particle number densities through the smoothed density field, \(\phi(\mathbf{r})\). With \(P\) being the cloud-in-cell (CIC) window function, the density and filtered densities are computed as

\[\phi(\mathbf{r}) = \sum_i P(\mathbf{r}-\mathbf{r}_i),\]

and

\[\tilde\phi(\mathbf{r}) = \int\mathrm{x}\mathbf{r}\, \phi(\mathbf{x})H(\mathbf{r}-\mathbf{x}),\]

where \(H\) is the grid-independent filtering function. The external potential is computed in reciprocal space as

\[V_\text{ext} = \mathrm{FFT}^{-1}\left[ \mathrm{FFT} \left( \frac{\partial w}{\partial \tilde\phi} \right) \mathrm{FFT}(H) \right],\]

where \(w\) is the interaction energy functional. Differentiating \(V_\text{ext}\) is done by simply applying \(i\mathbf{k}\) in Fourier space, and the resulting forces are back-transformed to direct space and interpolated to particle positions by

\[\mathbf{F}_j = -\sum_{j_k} \nabla V_{j_k} P(\mathbf{r}_{j_k} - \mathbf{r}_j)h^3,\]

where \(j_k\) are the neighbouring vertices of particle \(j\) at position \(\mathbf{r}_j\), and \(h^3\) is the volume of each grid voxel.

Parameters:

philist[pmesh.pm.RealField]

Pmesh RealField objects containing discretized particle number density values on the computational grid; one for each particle type. Pre-allocated, but empty; any values in this field are discarded. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_laplacianlist[pmesh.pm.RealField], (M, 3)

Like phi, but containing the laplacian of particle number densities.

phi_transferlist[pmesh.pm.ComplexField], (3,)

Like phi_fourier, used as an intermediary to perform FFT operations to obtain the gradient or laplacian of particle number densities.

layoutslist[pmesh.domain.Layout]

Pmesh communication layout objects for domain decompositions of each particle type. Used as blueprint by pmesh.pm.readout for exchange of particle information across MPI ranks as necessary.

force_meshlist[pmesh.pm.RealField]

Pmesh RealField objects containing discretized particle-field force density values on the computational grid; D fields in D dimensions for each particle type. Pre-allocated, but empty; any values in this field are discarded. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

hamiltonianHamiltonian

Particle-field interaction energy handler object. Defines the grid-independent filtering function, \(H\).

pmpmesh.pm.ParticleMesh

Pmesh ParticleMesh object interfacing to the CIC window function and the PFFT discrete Fourier transform library.

positions(N,D) numpy.ndarray

Array of positions for N particles in D dimensions. Local for each MPI rank.

types(N,) numpy.ndarray

Array of type indices for each of N particles. Local for each MPI rank.

configConfig

Configuration object.

v_extlist[pmesh.pm.RealField]

Pmesh RealField objects containing discretized particle-field external potential values on the computational grid; one for each particle type. Pre-allocated, but empty; any values in this field are discarded Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_fourierlist[pmesh.pm.ComplexField]

Pmesh ComplexField objects containing discretized particle number density values in reciprocal space on the computational grid; one for each particle type. Pre-allocated, but empty; any values in this field are discarded Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

v_ext_fourierlist[pmesh.pm.ComplexField]

Pmesh ComplesField objects containing discretized particle-field external potential values in reciprocal space on the computational grid; D+1 fields in D dimensions for each particle type. D copies are made after calculation for later use in force calculation, because the force transfer function application differentiates the field in-place, ruining the contents for differentiation in the remaining D-1 spatial directions. Pre-allocated, but empty; any values in this field are discarded. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

m: list[float], (M,)

pmesh.pm.ParticleMesh parameter for mass of particles in simulation unit. Defaults to 1.0 for all particle types.

compute_potentialbool, optional

If True, a D+1-th copy of the Fourier transformed external potential field is made to be used later in particle-field energy calculation. If False, only D copies are made.

See also

compute_field_and_kinetic_energy

Compute the particle-field energy after the external potential is calculated.

hymd.field.update_field_force_q(charges, phi_q, phi_q_fourier, psi, psi_fourier, elec_field_fourier, elec_field, elec_forces, layout_q, hamiltonian, pm, positions, config)

Calculate the electrostatic particle-field forces on the grid

Computes the electrostatic potential \(\Psi\) from particle charges through the smoothed charge density \(\tilde\rho\). With \(P\) being the cloud-in-cell (CIC) window function, the charge density and filtered charge densities are computed as

\[\rho(\mathbf{r}) = \sum_i q_i P(\mathbf{r}-\mathbf{r}_i),\]

and

\[\tilde\rho(\mathbf{r}) = \int\mathrm{x}\mathbf{r}\, \rho(\mathbf{x})H(\mathbf{r}-\mathbf{x}),\]

where \(H\) is the grid-independent filtering function. The electrostatic potential is computed in reciprocal space as

\[\Phi = \mathrm{FFT}^{-1}\left[ \frac{4\pi k_e}{\varepsilon \vert\mathbf{k}\vert^2} \mathrm{FFT}(\rho)\mathrm{FFT}(H) \right],\]

with the electric field

\[\mathbf{E} = \mathrm{FFT}^{-1}\left[ -i\mathbf{k}\,\mathrm{FFT}(\Psi) \right].\]

In the following, let \(\mathbf{F}_j\) denote the electrostatic force acting on the particle with index \(j\) and position \(\mathbf{r}_j\). The interpolated electrostatic force is

\[\mathbf{F}_j = \sum_k q_j\mathbf{E}_{j_k} P(\mathbf{r}_{j_k}-\mathbf{r}_j)h^3,\]

where \(\mathbf{E}_{j_k}\) is the discretized electric field at grid vertex \(j_k\), \(\mathbf{r}_{j_k}\) is the position of the grid vertex with index \(j_k\), and \(h^3\) is the volume of each grid voxel. The sum is taken over all closest neighbour vertices, \(j_k\).

Parameters:

charges(N,) numpy.ndarray

Array of particle charge values for N particles. Local for each MPI rank.

phi_qpmesh.pm.RealField

Pmesh RealField object for storing calculated discretized charge density density values on the computational grid. Pre-allocated, but empty; any values in this field are discarded. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_q_fourierpmesh.pm.ComplexField

Pmesh ComplexField object for storing calculated discretized Fourier transformed charge density values in reciprocal space on the computational grid. Pre-allocated, but empty; any values in this field are discarded. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

elec_field_fourierpmesh.pm.ComplexField

Pmesh ComplexField object for storing calculated discretized electric field values in reciprocal space on the computational grid. Pre-allocated, but empty; any values in this field are discarded. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

elec_fieldpmesh.pm.RealField

Pmesh RealField object for storing calculated discretized electric field values on the computational grid. Pre-allocated, but empty; any values in this field are discarded. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

elec_forces(N,D) numpy.ndarray

Array of electrostatic forces on N particles in D dimensions.

layout_qpmesh.domain.Layout

Pmesh communication layout object for domain decomposition of the full system. Used as blueprint by pmesh.pm.paint and pmesh.pm.readout for exchange of particle information across MPI ranks as necessary.

hamiltonianHamiltonian

Particle-field interaction energy handler object. Defines the grid-independent filtering function, \(H\).

pmpmesh.pm.ParticleMesh

Pmesh ParticleMesh object interfacing to the CIC window function and the PFFT discrete Fourier transform library.

positions(N,D) numpy.ndarray

Array of positions for N particles in D dimensions. Local for each MPI rank.

confighymd.input_parser.Config

Configuration object.

hymd.field.update_field_force_q_GPE(conv_fun, phi, types, charges, phi_q, phi_q_fourier, phi_eps, phi_eps_fourier, phi_eta, phi_eta_fourier, phi_pol_prev, phi_pol, elec_field, elec_forces, elec_field_contrib, psi, Vbar_elec, Vbar_elec_fourier, force_mesh_elec, force_mesh_elec_fourier, hamiltonian, layout_q, layouts, pm, positions, config, comm=<mpi4py.MPI.Intracomm object>)

Calculate the electrostatic particle-field forces on the grid, arising from a general Poisson equation, i.e. anisotropic permittivity/dielectric. The function is called when tomli input config.coulombtype = “PIC_Spectral_GPE.”

Computes the electrostatic potential \(\Psi\) from particle charges through the smoothed charge density \(\tilde\rho\). With \(P\) being the cloud-in-cell (CIC) window function, the charge density and filtered charge densities are computed as

\[\rho(\mathbf{r}) = \sum_i q_i P(\mathbf{r}-\mathbf{r}_i),\]

and

\[\tilde\rho(\mathbf{r}) = \int\mathrm{x}\mathbf{r}\, \rho(\mathbf{x})H(\mathbf{r}-\mathbf{x}),\]

where \(H\) is the grid-independent filtering function. The electrostatic potential for a variable dielectric does not have an analytical expression, and is computed in reciprocal through an iterative method.

The GPE states that

\[\nabla \cdot \left(\epsilon(\mathbf{r}) \nabla{\mathbf{\psi(r)}}\right) = -\rho({\mathbf{r}}).\]

where \(\epsilon(\mathbf{r})\) is the relative dielectric function.

Parameters:

conv_funConvergence function.

Returns a scalar. Depends on MPI allreduce for similar convergence across MPI ranks.

philist[pmesh.pm.RealField]

Pmesh RealField objects containing discretized particle number density values on the computational grid; one for each particle type. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

types(N,) numpy.ndarray

Array of type indices for each of N particles. Local for each MPI rank.

charges(N,) numpy.ndarray

Array of particle charge values for N particles. Local for each MPI rank.

phi_qpmesh.pm.RealField

Pmesh RealField object for storing calculated discretized charge density density values on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_q_fourierpmesh.pm.ComplexField

Pmesh ComplexField object for storing calculated discretized Fourier transformed charge density values in reciprocal space on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_epspmesh.pm.RealField

Pmesh RealField object for storing calculated discretized relative dielectric values on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_eps_fourierpmesh.pm.ComplexField

Pmesh ComplexField object for storing calculated discretized Fourier transformed relative dielectric values in reciprocal space on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_etapmesh.pm.RealField

Pmesh RealField object for storing calculated discretized gradients of the relative dielectric values on the computational grid. Pre-allocated,but empty. Changed in-place.Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_eta_fourierpmesh.pm.ComplexField

Pmesh ComplexField object for storing the calculated discretized Fourier transformed gradient relative dielectric values in reciprocal space on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_pol_prevpmesh.pm.RealField

Pmesh RealField object for storing calculated discretized polarization charge values on the computational grid. Parameter in the iterative method.Pre-allocated,but empty. Changedin-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

phi_polpmesh.pm.RealField

Pmesh RealField object for storing calculated discretized polarization charges on the computational grid. Parameter in the iterative method, updating the next quess in solving for the electrostatic potential. Pre-allocated,but empty. Changed in-place.Local for each MPI rank– the full computational grid is represented by the collective fields of all MPI ranks.

elec_fieldpmesh.pm.RealField

Pmesh RealField object for storing calculated discretized electric field values on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

elec_forces(N,D) numpy.ndarray

Array of electrostatic forces on N particles in D dimensions.

elec_field_contribpmesh.pm.RealField

Pmesh RealField object for storing \(|\mathbf{E(r)}|^2/\phi_{0}\) on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank– the full computational grid is represented by the collective fields of all MPI ranks.

psipmesh.pm.RealField

Pmesh RealField object for storing electrostatic potential on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank– the full computational grid is represented by the collective fields of all MPI ranks.

Vbar_elecmesh.pm.RealField

Pmesh RealField object for storing functional derivatives of :math:`|w({ phi })_{elec}`on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank– the full computational grid is represented by the collective fields of

all MPI ranks.

Vbar_elec_fourierpmesh.pm.ComplexField

Pmesh ComplexField object for storing the calculated functional derivatives of \(\|w(\{ \phi \})_{elec}\) in reciprocal space on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

force_mesh_elecpmesh.pm.RealField

Pmesh RealField object for storing electrostatic force values on the computational grid. Pre-allocated, but empty. Changed in-place. Local for each MPI rank– the full computational grid is represented by the collective fields of all MPI ranks.

force_mesh_elec_fourierpmesh.pm.ComplexField

Pmesh ComplexField object for storing the calculated electrostatic force values in reciprocal space on the computational grid. Local for each MPI rank–the full computational grid is represented by the collective fields of all MPI ranks.

hamiltonianHamiltonian

Particle-field interaction energy handler object. Defines the grid-independent filtering function, \(H\).

layout_qpmesh.domain.Layout

Pmesh communication layout object for domain decomposition of the full system. Used as blueprint by pmesh.pm.paint and pmesh.pm.readout for exchange of particle information across MPI ranks as necessary.

layouts: list[pmesh.domain.Layout]

Pmesh communication layout objects for domain decompositions of each particle type. Used as blueprint by pmesh.pm.readout for exchange of particle information across MPI ranks as necessary.

pmpmesh.pm.ParticleMesh

Pmesh ParticleMesh object interfacing to the CIC window function and the PFFT discrete Fourier transform library.

positions(N,D) numpy.ndarray

Array of positions for N particles in D dimensions. Local for each MPI rank.

confighymd.input_parser.Config

Configuration object.

comm: mpi4py.Comm

MPI communicator to use for rank commuication.

See also

compute_field_energy_q_GPE

Compute the electrostatic energy after electrosatic force is calculated for a variable (anisotropic) dielectric general Poisson equation.

4.8. File input/output

Handle file input/output in parllel HDF5 fashion

class hymd.file_io.OutDataset(dest_directory, config, double_out=False, disable_mpio=False, comm=<mpi4py.MPI.Intracomm object>)

HDF5 dataset handler for file output

close_file(comm=<mpi4py.MPI.Intracomm object>)

Closes the HDF5 output file

Parameters:

commmpi4py.Comm

MPI communicator to use for rank communication.

flush()

Flushes output buffers, forcing file writes

is_open(comm=<mpi4py.MPI.Intracomm object>)

Check if HDF5 output file is open

Parameters:

commmpi4py.Comm

MPI communicator to use for rank communication.

hymd.file_io.distribute_input(in_file, rank, size, n_particles, max_molecule_size=201, comm=<mpi4py.MPI.Intracomm object>)

Assign global arrays onto MPI ranks, attempting load balancing

Distributes approximately equal numbers of particles (workload) onto each independent MPI rank, while respecting the requirement that any molecule must be fully contained on a single MPI rank only (no splitting molecules across multiple CPUs).

Parameters:

in_fileh5py.File

HDF5 input file object.

rankint

Local rank number for this MPI rank.

sizeint

Global size of the MPI communicator (number of total CPUs).

n_particlesint

Total number of particles.

max_molecule_sizeint, optional

Maximum size of any molecule present in the system. Used to initially guess where the MPI rank boundaries (start/end indices) in the global arrays should be placed. If molecules of size >max_molecule_size exist in the simulation system, HyMD might work as expected. Or it might fail spectacularly.

commmpi4py.Comm

MPI communicator to use for rank commuication.

Returns:

rank_range

Starting and ending indices in the global arrays for each MPI rank.

hymd.file_io.setup_time_dependent_element(name, parent_group, n_frames, shape, dtype, units=None)

Helper function for setting up time-dependent HDF5 group datasets

All output groups must adhere to the H5MD standard, meaning a structure of

┗━ group particle group (e.g. all) or observables group

┗━ group time-dependent data

┣━ dataset step shape=(n_frames,)

┣━ dataset time shape=(n_frames,)

┗━ dataset value shape=(n_frames, *)

is necessary.

References

H5MD specification :

https://www.nongnu.org/h5md/h5md.html

hymd.file_io.store_data(h5md, step, frame, indices, positions, velocities, forces, box_size, temperature, pressure, kinetic_energy, bond2_energy, bond3_energy, bond4_energy, field_energy, field_q_energy, plumed_bias, time_step, config, velocity_out=False, force_out=False, charge_out=False, plumed_out=False, dump_per_particle=False, comm=<mpi4py.MPI.Intracomm object>)

Writes time-step data to HDF5 output file

Handles all quantities which change during simulation, as opposed to static quanitities (see store_static).

Parameters:

h5mdOutDataset

HDF5 dataset handler.

stepint

Step number.

frameint

Output frame number (step // n_print).

indices(N,) numpy.ndarray

Array of indices for N particles.

positions(N,) numpy.ndarray

Array of positions for N particles in D dimensions.

velocities(N,) numpy.ndarray

Array of velocities for N particles in D dimensions.

forces(N,) numpy.ndarray

Array of forces for N particles in D dimensions.

box_size(D,) numpy.ndarray

Array containing the simulation box size.

temperaturefloat

Calculated instantaneous temperature.

kinetic_energyfloat

Calculated instantaneous kinetic energy.

bond2_energyfloat

Calculated instantaneous harmonic two-particle bond energy.

bond3_energyfloat

Calculated instantaneous harmonic angular three-particle bond energy.

bond4_energyfloat

Calculated instantaneous dihedral four-particle torsion energy.

field_energyfloat

Calculated instantaneous particle-field energy.

field_q_energyfloat

Calculated instantaneous electrostatic energy.

plumed_biasfloat

PLUMED instantaneous bias energy.

time_stepfloat

Value of the time step.

configConfig

Configuration object.

velocity_outbool, optional

If True, velocities are written to output HDF5 file.

force_outbool, optional

If True, forces are written to output HDF5 file.

charge_outbool, optional

If True, electrostatic energies are written to the output HDF5 file.

plumed_outbool, optional

If True, PLUMED bias is written to the output HDF5 file.

dump_per_particlebool, optional

If True, all quantities are written per particle.

commmpi4py.Comm

MPI communicator to use for rank commuication.

See also

store_static

Outputs all static time-independent quantities to the HDF5 output file

hymd.file_io.store_static(h5md, rank_range, names, types, indices, config, bonds_2_atom1, bonds_2_atom2, molecules=None, velocity_out=False, force_out=False, charges=False, dielectrics=False, plumed_out=False, comm=<mpi4py.MPI.Intracomm object>)

Outputs all static time-independent quantities to the HDF5 output file

Parameters:

h5mdOutDataset

HDF5 dataset handler.

rank_rangelist[int]

Start and end indices for global arrays for each MPI rank.

names(N,) numpy.ndarray

Array of names for N particles.

types(N,) numpy.ndarray

Array of type indices for N particles.

indices(N,) numpy.ndarray

Array of indices for N particles.

configConfig

Configuration object.

bonds_2_atom1(B,) numpy.ndarray

Array of indices of the first particle for B total two-particle bonds.

bonds_2_atom2(B,) numpy.ndarray

Array of indices of the second particle for B total two-particle bonds.

molecules(N,) numpy.ndarray, optional

Array of integer molecule affiliation for each of N particles. Global (across all MPI ranks) or local (local indices on this MPI rank only) may be used, both, without affecting the result.

velocity_outbool, optional

If True, velocities are written to output HDF5 file.

force_outbool, optional

If True, forces are written to output HDF5 file.

charges(N,) numpy.ndarray

Array of particle charge values for N particles.

dielectrics(N,) numpy.ndarray

Array of particle relative dielectric values for N particles.

plumed_outbool, optional

If True, PLUMED bias is written to output HDF5 file.

commmpi4py.Comm

MPI communicator to use for rank commuication.

See also

prepare_bonds

Constructs two-, three-, and four-particle bonds from toplogy input file and bond configuration information.

distribute_input

Distributes input arrays onto MPI ranks, attempting load balancing.