Development / API reference

This section documents PIMMS’ internal Python API, generated automatically from the NumPy-style docstrings in the source. It is aimed at developers extending PIMMS or scripting against it; users configuring simulations want the keyword reference instead.

The package is organised around a few core objects: a Lattice holds the grid and the Chain objects; a Hamiltonian evaluates the energy; a MoveObject implements the Monte Carlo moves; an AcceptanceCalculator handles move selection and the Metropolis criterion; and the Simulation ties them together and drives the run, configured by the KeyFileParser.

Simulation engine

class pimms.simulation.Simulation(keyword_lookup)

Bases: object

The Simulation object: this is the master object from which all simulation-related events happen

ANAFUNCT_acceptance(step)

Run acceptance-criterion analysis.

Writes the current move acceptance statistics (held by the AcceptanceCalculator, self.ACC) to disk on each call.

Parameters:

step (int) – Current simulation step number, written alongside the statistics.

Return type:

None

ANAFUNCT_cluster_analysis(step)

Run cluster analysis.

Computes both the contact (short-range) and long-range cluster distributions for the current configuration, corrects cluster positions into a single periodic image, and derives polymeric properties, gross size/shape properties (volume, surface area, density) and radial density profiles for the size-thresholded clusters. All results are written to disk on each call, making this routine I/O heavy.

Parameters:

step (int) – Current simulation step number, written alongside the cluster data.

Return type:

None

ANAFUNCT_custom_stubb(step)

No-op custom-analysis stub.

Used as the ANA_CUSTOM analysis routine when no custom analysis module is side-loaded via the keyfile. Does nothing.

Parameters:

step (int) – Current simulation step number (ignored).

Return type:

None

ANAFUNCT_distance_map(step)

Run distance map analysis.

Updates the running distance-map counters associated with each chain on the lattice. No data is written to disk on each call; the accumulated map is written at the end of the simulation.

Parameters:

step (int) – Current simulation step number (accepted for interface uniformity; not used directly).

Return type:

None

ANAFUNCT_end_to_end(step)

Run end-to-end distance analysis.

Computes the end-to-end distance for every chain and writes the list to disk on each call, making this routine I/O heavy.

Parameters:

step (int) – Current simulation step number, written alongside the data.

Return type:

None

ANAFUNCT_internal_scaling(step)

Run internal scaling analysis.

Updates the running internal-scaling counters (both normal and squared) associated with each chain on the lattice. No data is written to disk on each call; the accumulated averages are written at the end of the simulation.

Parameters:

step (int) – Current simulation step number (accepted for interface uniformity; not used directly).

Return type:

None

ANAFUNCT_polymeric_properties(step)

Run polymeric-properties (radius of gyration and asphericity) analysis.

Computes the radius of gyration and asphericity for every chain and writes both lists to disk on each call, making this routine I/O heavy.

Parameters:

step (int) – Current simulation step number, written alongside the data.

Return type:

None

ANAFUNCT_save_restart(step)

Write a restart file capturing the current system state.

Builds a restart.RestartObject from the current lattice (recording the hardwall status), stores the freshly evaluated total energy in it, and writes the restart file to disk. The resulting file can be used to re-initialize the system in a later simulation.

Parameters:

step (int) – Current simulation step number, used only for the status message.

Return type:

None

auxillary_chain_update(old_energy)

Advance (and possibly finalize) an in-progress system-wide TSMMC move.

Performs all the busywork associated with the system-wide temperature switch Metropolis Monte Carlo (TSMMC) move. When the temperature sweep held by self.TSMMC_coordinator is not yet complete this simply checks the auxiliary chain in (updating the acceptance object only if the temperature changed). When the sweep is complete it applies the accept/reject decision: on rejection the lattice is restored from the coordinator’s backup, on acceptance the (already-updated) lattice is kept.

Whether this logic should live here, inside the TSMMC object, or inside the MOVER object is undecided; for now it lives here on the logic that it makes global changes to the simulation system and so belongs with the Simulation object.

Parameters:

old_energy (int or float) – Current total system energy, used both for the acceptance test and for reporting the energy change of a completed move.

Returns:

A two-place tuple (move_complete, move_accepted). move_accepted is always False while the move has not yet completed.

Return type:

tuple of (bool, bool)

build_R2R_distance_distribution_analysis(R2R_info)

This function returns a function which is initialized by the variables pass in by the R2R_info - in essence generating a closure.

Basically, if you’re not familiar with functional programming, this creates a new function where the $R2R_info variable inside the ANAFUNCT_R2R_distance function is set by the build_R2R_distance_distribution_analysis function.

This function (ANAFUNCTION_R2R_distance) is then returned, and next time its called the R2R_info variable IN THE FUNCTION BEING CALLED is already initialzed.

Parameters:

R2R_info (list) – List of residue-index pairs (i, j) for which the residue-residue distance distribution should be computed and written. May be empty.

Returns:

A closure ANAFUNCT_R2R_distance(step) that, when called, computes the requested residue-residue distances across all chains for the given step and writes them to disk.

Return type:

callable

end_of_simulation_analysis()

Final analysis routines run at the end of the simulation. For all analysis where a final average value makes sense this is going to be where the code to calculate and save that output is written.

Currently computes and writes the chain-averaged internal scaling, scaling exponents (nu, R0) and distance maps, handling both single-chain-type and multicomponent (per-chain-type) systems.

Return type:

None

quench_update(i, old_energy)

Apply a temperature-quench update on quench steps.

Helper that runs a quench update if a temperature quench run is being performed. On steps that are multiples of self.QUENCH_FREQ it either reports (and disables the quench) when the target temperature has been reached, or advances the temperature one QUENCH_STEPSIZE toward the target (negating the step for heating quenches). When TSMMC is in use the TSMMC_coordinator is rebuilt at the new temperature, and the quench event is written to QUENCH.dat. Updates all relevant simulation variables in place.

Parameters:

• i (int) – Current simulation step number.

• old_energy (int or float) – Current total system energy, written to the quench output file when a temperature change occurs.

Return type:

None

Raises:

SimulationException – If self.QUENCH_FREQ is not a positive integer.

rigid_cluster_move(new_chain_positions, old_chain_positions)

Function which implements optimized energy calculations for rigid body cluster moves.

---

NOTE that unlike the single chain moves we actually determine the set of moved pairs inside this function. There’s a reason for this! So, when making a rigid cluster move we first determine the positions of all the chains in the cluster we’re moving. This provides us with a useful set of prior information because we KNOW that in that clusters’ original position the ONLY short range interactions we care about are between lattice sites occupied by the cluster and lattice sites occupied by the solvent. We know this because any sites between a cluster-component and a NON solvent site would be an intra-cluster pair, and given rigid cluster movements cannot be changing intra-cluster interactions the change in energy associated with intracluster sites must be zero. However, the key computational cost here is evaluating how the long range interactions contribute to the cost of moving the cluster.

Because we know all the relative interactions WITHIN the cluster must be held fixed (both short and long-range) the only interactions we care about are between the cluster interface. If we have NO long range interactions we can actually just perform the move without worrying about the energy because - by definition - we cannot be moving a cluster into direct contact with another solute molecule so the cluster-system interface is purely solute-solvent before and after - i.e. no change in energy. If we do have long range interactions (as would be usual) we have to compute their influence.

1. Determine all the long-range interactions goin’ on

2. Just TRANSLATE/ROTATE those positions to get the interfacial pairs in the clusters’ new position

However, to do this we need the cluster back in its original position - hence why we have to move the lattice BACK to its original position before we determine the interfacial residues

---

Parameters:

• new_chain_positions (dict) – Mapping of chainID to the list of new positions for that chain. The full set of new positions for the chains making up the cluster.

• old_chain_positions (dict) – Mapping of chainID to the list of original positions for that chain. The full set of old positions for the chains making up the cluster.

Returns:

The change in total system energy produced by the rigid cluster move. Returns 0.0 immediately for the energy-neutral case where the Hamiltonian has no long-range interactions (the move is still committed to the grids in that case).

Return type:

float

Raises:

Exception – If new_chain_positions and old_chain_positions do not describe the same set of chainIDs.

rigid_cluster_revert(new_chain_positions, old_chain_positions)

Revert a rejected rigid cluster move.

Restores the system back to its pre-move state after a rigid cluster move is rejected. Every chain in the cluster is deleted from its new positions (on both the grid and the type_grid) and then re-inserted at its original positions, and each chain object’s ordered positions are reset.

Parameters:

• new_chain_positions (dict) – Mapping of chainID to the list of new (rejected) positions for that chain.

• old_chain_positions (dict) – Mapping of chainID to the list of original positions to restore for that chain.

Return type:

None

run_all_analysis(step)

Master analysis function - cycles over each type of analysis to assess if that analysis should be performed this step (or not), and launching the analysis if it should be done.

Updates various analysis state information (e.g. internal scaling distances associated with each Chain etc.) but does not change any lattice positions or anything like that

Analysis is skipped entirely while the simulation is still in equilibration (step < self.equilibration). Non-default-frequency analysis routines run on their own per-routine frequencies, while default-frequency routines run together every self.anafreq steps.

Parameters:

step (int) – Current simulation step number, used to decide which analysis routines (if any) fire this step.

Return type:

None

run_simulation()

Execute the full Monte Carlo simulation workflow.

This method drives a complete production run from the current Simulation state. It performs one-time startup tasks, iterates over self.n_steps Monte Carlo steps, dispatches move proposals, applies Metropolis acceptance/rejection logic, runs analysis and trajectory/energy I/O at configured frequencies, and performs final output and cleanup.

High-level flow

1. Record the global start time and print startup status.

2. Compute initial system energy with the full Hamiltonian.

3. Initialize quench and trajectory output files when enabled.

4. Run startup analysis/file initialization routines.

5. Enter the main simulation loop and repeat until n_steps is

   reached:

   • Skip move proposals entirely if all chains are frozen.

   • If not in an auxiliary TSMMC chain:

     • apply quench updates (if QUENCH_RUN),

     • handle box-resize equilibration logic (if resize_eq),

     • perform non-analysis simulation I/O,

     • run scheduled analysis callbacks.

   • If inside an auxiliary TSMMC chain, update/complete that chain

     and only advance the global step counter once the TSMMC cycle is complete.

   • Choose a random movable chain and sample a move type via

     AcceptanceCalculator.move_selector.

   • Execute the selected move implementation (single-chain,

     cluster, TSMMC, system-shake, etc.).

   • For standard move families, compute energy deltas and apply

     Boltzmann acceptance; on rejection, revert lattice state.

   • Update move statistics used for post-hoc diagnostics.

6. After the loop, optionally flush an in-memory trajectory (when

   SAVE_AT_END is active), run final analysis, and always write a final restart snapshot.

Side effects

• Mutates simulation state in place, including lattice coordinates,

  chain positions, energies, counters, and acceptance statistics.

• Writes multiple output artifacts (trajectory, energy, quench,

  analysis files, restart file), depending on runtime options.

• Emits progress and diagnostic messages to stdout/loggers.

Notes

• Invalid move-selection codes raise SimulationException.

• Energy consistency checks may raise

  SimulationEnergyException via simulation_IO.

• The method is intentionally monolithic because step ordering is

  coupled to detailed-balance constraints and output semantics.

rtype:

None

setup_analysis(keyword_lookup)

This function constructs two dictionaries. Each key-value pair in the dictionary is a function-frequency pair, where the function is an analysis routine of the format FXC(step) and the frequency is the frequency with which that analysis is performed.

The two lists correspond to analysis which occurs with the same frequency as the general analysis and then the analysis which occurs at a frequency different to the general analysis.

This is a bit of work at the start, but allows us to run bespoke, custom-frequency analysis in a very simple way during the simulation.

keyword_lookup provides all the info needed.

Parameters:

keyword_lookup (dict) – The controlled-vocabulary keyword dictionary. Provides each analysis keyword’s frequency, the residue pairs for R2R analysis (ANA_RESIDUE_PAIRS), the default frequency (ANALYSIS_FREQ), and an optional side-loaded custom ANALYSIS_MODULE.

Returns:

(non_default_freq_analysis, default_freq_analysis). Each is a mapping from an analysis function (a bound method or closure taking a single step argument) to the integer frequency at which it should run. The first holds routines whose frequency differs from the default analysis frequency; the second holds those that match it.

Return type:

tuple of (dict, dict)

Raises:

SimulationException – If the internal failsafe consistency checks on the analysis keyword tables fail (indicates a software bug).

simulation_IO(i, old_energy)

Helper function to run simulation IO (O) for various different things. Additional output should be added her as an if statement comparing against the appropriate keyword frequency.

NOTE: All analysis is dealt with seperatly and shouldn’t be added here - this is for non-analysis IO (i.e. status IO).

This includes

1. Printing status of the simulation

2. Writing out trajectory information

3. Writing out energy information

4. Performing global energy comparison

Parameters:

• i (int) – Current step that the simulation is on

• old_energy (int or float) – Current (locally-tracked) total system energy. Used for status printing, written to the energy file, and compared against the from-scratch Hamiltonian recalculation during the periodic energy consistency check.

Return type:

None

Raises:

SimulationEnergyException – If, on an energy-comparison step, the locally-tracked energy differs from the fully recalculated energy (a configuration snapshot is written to CONFIG_AT_ENERGY_FAIL.pdb/.xtc before raising).

single_chain_move(move_event, chainID)

Compute the energy change of a single-chain move.

Implements the optimized local energy calculation for moves that perturb a single chain. The method evaluates the short-range, long-range, super long-range and angle contributions for the chain in both its original and moved positions (operating only on the local interaction envelope rather than the whole lattice) and commits the chain to its new position on the grid, type_grid and chain object. The returned value is the resulting change in total system energy, which the caller passes to the Metropolis acceptance test.

Parameters:

• move_event (MoveEvent) – Object containing all the move details (original/moved positions, moved indices, move type, etc.).

• chainID (int) – ID of the single chain being moved. Each chain has a unique ID starting at 1 and increasing.

Returns:

The change in total system energy (local_dif) produced by the move, including short-range, long-range, super long-range and angle terms.

Return type:

float

single_chain_revert(move_event, chainID)

Revert a rejected single-chain move.

Restores the system back to its pre-move state after a single-chain move is rejected by the Metropolis criterion. The chain is removed from its moved positions and placed back at its original positions on the grid, the chain object’s ordered positions are reset, and the type_grid is updated back.

Parameters:

• move_event (MoveEvent) – Object containing the move details (moved/original positions and chain positions, and moved indices) used to undo the move.

• chainID (int) – ID of the single chain whose move is being reverted.

Return type:

None

startup_analysis()

Function for including all the analysis activity which should be run BEFORE the simulation starts. This includes wiping files we’re going to progressively append to during the simulation, but one could imagine in the future additional code might be included here for functional purposes.

For multicomponent systems (more than one chain type) it additionally wipes per-chain-type cluster output files.

Return type:

None

update_dimensions(step, old_energy)

Resize the lattice box at the end of resize-equilibration.

Handles the box-resize-equilibration logic. On any step other than the end of equilibration this is a no-op that returns the energy unchanged. On the final equilibration step it checks whether any chain still straddles the periodic boundary:

• If one or more chains straddle the boundary, a warning is logged, the number of steps and the equilibration length are each extended by 100, and the offending chains are returned as a forced-move override list.

• Otherwise the lattice is rebuilt at the production dimensions via a restart.RestartObject (optionally applying EQ_OFFSET), output trajectory/PDB files are (re)initialized, the resize flag is cleared, hardwall is switched off if the production run uses PBC, and the energy is recomputed from scratch.

Parameters:

• step (int) – Current simulation step number, compared against self.equilibration.

• old_energy (int or float) – Current total system energy, returned unchanged except when the box is actually resized (in which case the energy is recomputed).

Returns:

(chain_selection_override, energy). chain_selection_override is an empty list except when chains still straddle the boundary, in which case it is the list of offending chainIDs that must be forced to move. energy is the (possibly recomputed) total system energy.

Return type:

tuple

Monte Carlo moves

class pimms.acceptance.AcceptanceCalculator(temp, keyword_lookup)

Bases: object

The AcceptanceCalculator object is a object associated with each simulation (though if Hamiltonian switch MC is implemented later then you’d probably want one AcceptanceCalculator per Hamiltonian)

The object performs the following functions

1. Keeps track of the simulation temperature

2. Implements move acceptance/rejection (boltzmann_acceptance)

3. Record move attempts and success for post-hoc analysis (update_move_logs()

4. Perform the actual randomized move selection (move_selector())

alt_Markov_chain_update_move_logs(number_tried)

Update the counter tracking proposed moves made in alternative Markov chains.

These are the accept/reject events performed inside various sub-moves (e.g. the TSMMC temperature-excursion moves) that belong to an alternative Markov chain rather than the main one. This bookkeeping is for performance reasons only (i.e. correctly comparing the number of independent accept/reject events between different simulation strategies). If running inside an auxiliary TSMMC chain the parallel aux_chain_alt_Markov_chain_moves counter is updated instead.

Parameters:

number_tried (int) – The number of alternative-Markov-chain moves to add to the counter.

Return type:

None

boltzmann_acceptance(old_energy, new_energy)

Accept or reject a move based on dE and the Boltzmann acceptance criterion.

Downhill (or energy-neutral) moves are always accepted. Uphill moves are accepted with probability exp(-(new_energy - old_energy) * invtemp), compared against a uniform random draw.

Parameters:

• old_energy (float) – The system energy before the proposed move.

• new_energy (float) – The system energy after the proposed move.

Returns:

True if the move is accepted, False otherwise.

Return type:

bool

get_total_aux_chain_moves()

Calculate the total number of moves attempted in auxiliary Markov chains.

Sums every per-move auxiliary-chain attempt counter (aux_chain_move_count) together with the auxiliary alternative Markov-chain move counter (aux_chain_alt_Markov_chain_moves). Note this considers only the number of moves attempted, not accept/reject information.

Returns:

The total number of moves attempted across all auxiliary Markov chains.

Return type:

int

megastep_update_move_logs(selection, number_accepted, number_tried)

Bulk-update the move logs for a megamove that performed many sub-moves.

Like update_move_logs(), but increments the attempt and accepted counters by arbitrary amounts in a single call - used by the megamoves (crankshaft / slither / pull) which internally perform and accept/reject many sub-moves per call. move_count and accepted_count are object lists where selection defines which move is being updated (1 = crankshaft, 2 = chain translate, etc); the 0th element is ignored. If running inside an auxiliary TSMMC chain the parallel aux_chain_* counters are updated instead.

Parameters:

• selection (int) – The MoveType code of the move being logged.

• number_accepted (int) – The number of sub-moves that were accepted.

• number_tried (int) – The number of sub-moves that were attempted/proposed.

Return type:

None

Raises:

AcceptanceException – If selection is out of the valid 1..NUM_MOVES range.

move_selector(chain_length)

Based on the MOVESET defined the keyfile, returns a value between 1 and 8 which corresponds to a specific move type, as defined below. Selection of these numbers is defined based on the frequency specificed in the keyfile by the MOVE_* parameters

If the chain selected is of length 1 then the moveset defaults to either a chain translation (rotation makes no sense) OR a cluster rotation or translation. This behaviour is actually hardcoded - specifically because if you have a system with various chains then this allows optimal move behaviour for single-residue chains vs. multiresidue chains.

In the future it might be wise to allow each chain-type to have a move set defined with it.

CODE | MOVE ______________________________ 1 | CRANKSHAFT 2 | CHAIN TRANSLATE 3 | CHAIN ROTATE 4 | CHAIN PIVOT 5 | HEAD PIVOT 6 | SLITHER 7 | CLUSTER TRANSLATE 8 | CLUSTER ROTATE 9 | TSMMC CHAIN re-arrangement 10 | TSMMC MULTICHAIN re-arrangement 11 | PULL 12 | TSMMC SYSTEM re-arrangement 13 | JUMP AND RELAX 14 | VMMC

Parameters:

chain_length (int) – The length of the chain selected to be moved. If this is 1, moves that are meaningless for a single bead (chain rotate, chain pivot, head pivot, slither and pull; codes 3, 4, 5, 6, 11) are remapped to a crankshaft move (code 1).

Returns:

The selected MoveType code (an integer between 1 and 14). If the calculator is currently operating inside an auxiliary TSMMC chain (self.auxillary_chain is True), the nested TSMMC moves (codes 9, 10, 12) are remapped to a crankshaft move to avoid nesting TSMMC excursions.

Return type:

int

Raises:

AcceptanceException – If no valid move could be selected, which indicates a bug in how the move-selection thresholds were constructed.

total_attempted_moves()

Total number of individual accept/reject events attempted across ALL sub-loops so far (i.e. every MC move, not just outer master-loop steps).

Megamoves (crankshaft / slither / pull) contribute their full per-megamove substep counts, single-chain moves contribute one each, and the TSMMC temperature-excursion moves contribute their alternative-Markov-chain substeps. This is the same quantity written to TOTAL_MOVES.dat.

Returns:

The cumulative number of attempted MC moves.

Return type:

int

update_move_logs(selection, acceptance)

Update the log of which moves are attempted and which are accepted.

move_count and accepted_count are object lists where selection defines which move is being updated (1 = crankshaft, 2 = chain translate, etc). There is a 0th element, but it is just ignored. If the calculator is currently running inside an auxiliary TSMMC chain (self.auxillary_chain), the parallel aux_chain_* counters are updated instead.

Parameters:

• selection (int) – The MoveType code of the move being logged.

• acceptance (bool) – True if the move was accepted (increments the accepted counter), False otherwise (only the attempt counter is incremented).

Return type:

None

Raises:

AcceptanceException – If selection is out of the valid 1..NUM_MOVES range.

update_temperature(temp)

Update the acceptance object’s temperature dynamically.

Re-validates the new temperature and recomputes the cached inverse temperature (self.invtemp). Used by quenching and by the TSMMC temperature-excursion machinery.

Parameters:

temp (float) – The new temperature (must be > 0).

Return type:

None

Raises:

AcceptanceException – If temp is not strictly greater than zero.

Energy

class pimms.energy.EmptyHamiltonian

Bases: object

Dummy class which implements a Hamiltonian used when all interactions are turned off. Means we can define whatever stub-functionality here without adding code-rot to the true Hamiltonian class

convert_sequence_to_LR_integer_sequence(sequence)

Stub long-range residue-to-integer conversion for the non-interacting case.

Mirrors Hamiltonian.convert_sequence_to_LR_integer_sequence() but, because no residue undergoes long-range interactions, always returns an empty list.

Parameters:

sequence (list or str) – The (human readable) residue sequence to convert. Ignored.

Returns:

Always an empty list.

Return type:

list

convert_sequence_to_integer_sequence(sequence)

Stub residue-to-integer conversion for the non-interacting case.

Mirrors Hamiltonian.convert_sequence_to_integer_sequence() but, because there are no interactions, simply maps every residue to the same dummy integer code (1).

Parameters:

sequence (list or str) – The (human readable) residue sequence to convert.

Returns:

A list of 1 values with the same length as sequence.

Return type:

list of int

evaluate_angle_energy(x, y)

Stub angle-energy evaluation that always returns zero.

Parameters:

• x (object) – Placeholder argument (typically chain positions). Ignored.

• y (object) – Placeholder argument (typically an intcode sequence). Ignored.

Returns:

Always 0.0.

Return type:

float

evaluate_local_energy(x, y)

Stub short-range local-energy evaluation that always returns zero.

Parameters:

• x (object) – Placeholder argument (typically a lattice object). Ignored.

• y (object) – Placeholder argument (typically a list of pairs). Ignored.

Returns:

Always 0.0.

Return type:

float

evaluate_local_energy_LR(x, y)

Stub long-range local-energy evaluation that always returns zero.

Parameters:

• x (object) – Placeholder argument (typically a lattice object). Ignored.

• y (object) – Placeholder argument (typically a list of pairs). Ignored.

Returns:

Always 0.0.

Return type:

float

evaluate_total_energy(x)

Stub total-energy evaluation that always returns zero.

Parameters:

x (object) – Placeholder argument (typically a lattice object). Ignored.

Returns:

Always 0.0.

Return type:

float

get_indices_of_long_range_residues(sequence)

Stub lookup of long-range residue indices for the non-interacting case.

Mirrors Hamiltonian.get_indices_of_long_range_residues() but, because no residue undergoes long-range interactions, always returns an empty list.

Parameters:

sequence (list or str) – The (human readable) residue sequence to inspect. Ignored.

Returns:

Always an empty list.

Return type:

list

class pimms.energy.Hamiltonian(parameter_file, num_dimensions, non_interacting, angles_off, hardwall=False, temperature=False, reduced_printing=False)

Bases: object

Main Hamiltonian class for evaluating system energies

build_angle_interactions(angle_dict, num_dimensions, angles_off)

Function that constructs a lookup table that we use to assign ‘angle’ withstraints. The angle effects are really related to the 1_3 interaction, but can also be used

The resulting integer-typed lookup table is stored on the instance as self.angle_lookup. Its shape depends on dimensionality: in 3D it is (max_intcode + 1, 3, 3, 3, 3, 3, 3) and in 2D it is (max_intcode + 1, 3, 3, 3, 3). The first axis is indexed by the residue integer code and the remaining axes encode the relative offsets of the i-1 and i+1 neighbours. Each residue must have an associated angle definition (three penalties A1/A2/A3); a missing definition raises an exception. Non-integer penalties are rounded to the nearest integer because the lattice energies are integer valued. When angles are switched off (or when only solvent exists) the table is populated with zeros.

Parameters:

• angle_dict (dict or bool) – Mapping from residue name to its [A1, A2, A3] angle penalties, as parsed from the parameter file. When angles_off is True this argument is ignored (and is typically passed as False).

• num_dimensions (int) – The lattice dimensionality. Must be 2 or 3.

• angles_off (bool) – If True, all angle penalties are forced to zero and the residue names are taken from the interaction table rather than from angle_dict.

Returns:

The function does not return a value; it sets self.angle_lookup as a side effect.

Return type:

None

Raises:

• EnergyException – If num_dimensions is neither 2 nor 3.

• ParameterFileException – If a residue has through-space interactions defined but no corresponding angle energies.

build_interaction_table(non_interacting=False)

Carries out dynamic construction of a two AxA float matrix where indicies along X and Y axis correspond to residues defined in the parameter file.

• The residue interaction table (RIT) defines short-range interactions

• The LRRIT (long range interaction table) defines long-range interactions

The return values are as follows:

RIT - 2D numpy array of floats which describes the short range interactions

The matrix is indexed using integer codes, where each code maps to a specific residue type

LRRIT - 2D numpy array of floats which describes the long range interactions

The matrix is indexed using integer codes, where each code maps to a specific residue type (same mapping as the RIT). The LR_RIT and the RIT are the same size (which allows the same indexing codes to be used) but MANY residues will not engage in LR interactions, so those sites are set to np.NaN such that if a bug leads to these being used an error will occur.

MAPPING - a residue-to-code mapping allowing for a residue name to be mapped

to its integer code (same code for short range and long range). This is a dictionary, where the key is the residue name (string) and the value is the corresponding integer code

LR_MAPPING - mapping where ONLY LR residues are included. This is a subset of

the mapping defined in MAPPING but has the nice feature of only including LR residues (so also offers a simple way to determine if a given residue undergoes LR interactions or not).

SLRRIT - 2D numpy array of floats which describes the super long range

interactions (SLR). The matrix is indexed using integer codes, where each code maps to a specific residue type (same mapping as the RIT and the LR_RIT).

Parameters:

non_interacting (bool) – If True, every entry in the short-range, long-range and super-long-range tables is forced to zero (overriding the parameter file), and a warning is emitted per residue pair unless reduced printing is enabled. Default is False.

Returns:

A 5-tuple (RIT, MAPPING, LRRIT, LR_MAPPING, SLRRIT) as described above: the short-range interaction table, the residue-name-to-integer mapping, the long-range interaction table, the long-range-only residue-name-to-integer mapping, and the super-long-range interaction table.

Return type:

tuple

convert_sequence_to_LR_integer_sequence(sequence)

This takes an human name residue sequence (e.g. amino acid) and converts it into it’s local LONG RANGE integer-code sequence. There is a unique mapping between each amino acid and it’s integer LR residue code. This is dynamically constructed on each simulation run and should not be assumed to be the same between different simulations.

Also important, the LR_int code is not the same as the short range int_code. The LR_int code is defined by the self.LR_parameter_to_int_map, while the int code is defined by the self.parameter_to_int_map.

Also importantly, EVERY residue must have a parameter_to_int_map entry but not every needs as LR_parameter_to_int_map entry.

Basically, this happens because Cython can’t do lookups with strings/chars but can with INTs, so we convert a list of chars into a list of ints

Parameters:

sequence (list or str) – The (human readable) residue sequence to convert.

Returns:

The long-range integer-code sequence, one entry per residue. A value of -1 is used for any residue that does not participate in long-range interactions (i.e. is absent from self.LR_parameter_to_int_map).

Return type:

list of int

convert_sequence_to_integer_sequence(sequence)

This takes an human name residue sequence (e.g. amino acid) and converts it into it’s local integer-code sequence which is then evaluated by the HYPERLOOP code - basically this happens because CYTHON can’t do lookups with strings/chars but can with INTs, so we convert a list of chars into a list of ints and then go to TOWN on that badboy

Parameters:

sequence (list or str) – The (human readable) residue sequence to convert. Each element must be a key in self.parameter_to_int_map.

Returns:

The integer-code sequence, one integer per residue in sequence.

Return type:

list of int

Raises:

ParameterFileException – If a residue in sequence has no corresponding integer code in the parameter file mapping.

evaluate_angle_energy(chain_positions, intcode_sequence, dimensions)

Angle energies are determined by hyperloop functions that basically compare the angle vector and use a pre-computed lookup table to convert that angle into some energy penalty. It’s lightning fast, and residue specific!

Chains with fewer than 3 beads have no defined angle and contribute zero. For valid chains the work is dispatched to the 2D or 3D Cython angle hyperloop (evaluate_angle_energy_2D / evaluate_angle_energy_3D) using the pre-computed self.angle_lookup table.

Parameters:

• chain_positions (list) – A list of 2D or 3D bead positions (each a 2- or 3-element list of integer lattice coordinates) over which the angle energy is evaluated.

• intcode_sequence (list of int) – An equal-length list of the per-bead integer codes used to index the residue-specific angle lookup table.

• dimensions (list) – A list of the box dimensions; its length (2 or 3) determines the lattice dimensionality.

Returns:

The total angle penalty for the supplied chain. Returns 0.0 when the chain has fewer than 3 beads.

Return type:

float or int

Raises:

EnergyException – If the lattice dimensionality is neither 2 nor 3.

evaluate_local_energy(latticeObject, pairs_list)

This is really the main energy calculating function - it takes a latticeObject (which contains the comple information on what residues are where) and a pairs_list which defines the pairs of residues that the energy will be calculated over (i.e. defining the ‘locality’ of this operation - LOCAL does not here mean only short range!).

This particular method evaluates the SHORT-RANGE energy by dispatching to __evaluate_local_energy_shortrange using the short-range residue interaction table.

Parameters:

• latticeObject (Lattice) – The lattice object holding the type grid and box dimensions.

• pairs_list (numpy.ndarray or list) – The (non-redundant) set of position pairs over which the short-range energy is evaluated.

Returns:

The total short-range interaction energy over the supplied pairs. Returns 0 when pairs_list is empty.

Return type:

float or int

evaluate_local_energy_LR(latticeObject, pairs_list)

This is really the main energy calculating function - it takes a latticeObject (which contains the comple information on what residues are where) and a pairs_list which defines the pairs of residues that the energy will be calculated over (i.e. defining the ‘locality’ of this operation - LOCAL does not here mean only short range!).

This particular method evaluates the LONG-RANGE energy by dispatching to __evaluate_local_energy_non_shortrange using the long-range residue interaction table.

Parameters:

• latticeObject (Lattice) – The lattice object holding the type grid and box dimensions.

• pairs_list (numpy.ndarray or list) – The (non-redundant) set of position pairs over which the long-range energy is evaluated.

Returns:

The total long-range interaction energy over the supplied pairs. Returns 0 when pairs_list is empty.

Return type:

float or int

evaluate_local_energy_SLR(latticeObject, pairs_list)

This is really the main energy calculating function - it takes a latticeObject (which contains the comple information on what residues are where) and a pairs_list which defines the pairs of residues that the energy will be calculated over (i.e. defining the ‘locality’ of this operation - LOCAL does not here mean only short range!).

This particular method evaluates the SUPER-LONG-RANGE (SLR) energy by dispatching to __evaluate_local_energy_non_shortrange using the super-long-range residue interaction table.

Parameters:

• latticeObject (Lattice) – The lattice object holding the type grid and box dimensions.

• pairs_list (numpy.ndarray or list) – The (non-redundant) set of position pairs over which the super-long-range energy is evaluated.

Returns:

The total super-long-range interaction energy over the supplied pairs. Returns 0 when pairs_list is empty.

Return type:

float or int

evaluate_total_energy(latticeObject, id_to_typeMap=None)

Function which evaluates the total energy of the system.

This function constructs a redundant list of ALL pairwise interactions between each residue and every neighbour, and then uses the evaluate_local_energy to define the energy of the ‘local’ system defined by those pairs (where the local system happens to actually be the entire system)

Concretely, it walks every chain on the lattice, gathers the ordered bead positions and the per-bead long-range binary array, accumulates the angle energy chain-by-chain, builds the non-redundant short-range, long-range and super-long-range pair lists via lattice_utils.build_all_envelope_pairs and then evaluates each energy contribution.

Parameters:

• latticeObject (Lattice) – The lattice object holding all chains, the type grid and the box dimensions over which the total energy is computed.

• id_to_typeMap (dict, optional) – Currently unused. Retained for interface compatibility. Default is None.

Returns:

A 5-tuple (total, energy_local, energy_LR, energy_SLR, angle_energy) where total is the sum of the short-range (local), long-range, super-long-range and angle energy contributions, and the remaining elements are those individual contributions.

Return type:

tuple of float

get_indices_of_long_range_residues(sequence)

Return the sequence indices of residues that undergo long-range interactions.

Takes an amino acid sequence and returns a list with the index of positions which undergo long-range interactions. A residue is considered long-range if it appears as a key in self.human_readable_LR_interaction_table.

Parameters:

sequence (list or str) – The (human readable) residue sequence to inspect.

Returns:

The (zero-based) indices into sequence of residues that engage in long-range interactions.

Return type:

list of int

set_hardwall(value=True)

Toggle whether interactions use hardwall or periodic boundaries.

Internally the hardwall flag is stored as an integer (1 for hardwall, 0 for periodic boundary conditions) because it is passed through to the Cython energy kernels.

Parameters:

value (bool) – If True the Hamiltonian uses hardwall boundary conditions (self.hardwall = 1); if False it uses periodic boundary conditions (self.hardwall = 0). Default is True.

Return type:

None

Lattice, chains and geometry

pimms.lattice_utils.append_to_xtc_file(lattice, spacing, xtc_filename='traj.xtc', autocenter=False)

Low level function that adds a current lattice to and existing XTC file

Parameters:

• lattice (lattice.Lattice) – A lattice object

• spacing (float) – Lattice-to-realspace spacing in angstroms.

• xtc_filename (str) – Filename to read from and extend

• autocenter (bool) – Flag which, if set to True and there’s a single chain will center the protein in the box. This is useful for visualization purposes but does mean any translational diffusion will be lost. Default = False

Returns:

No return by the existing XTC file is extended by one frame

Return type:

None

pimms.lattice_utils.append_to_xtc_file_non_redundant(lattice, spacing, pdb_filename='START.pdb', xtc_filename='traj.xtc', autocenter=False, unwrap=False)

Low level function that adds a current lattice to an existing XTC file. This is different than ‘append_to_xtc_file’ in that it does not make a frame.pdb object each time it needs to save the frame. Uses the START.pdb file as the topology.

Parameters:

• lattice (lattice.Lattice) – A lattice object

• spacing (float) – Lattice-to-realspace spacing in angstroms.

• pdb_filename (str) – Topology filename to read from and extend

• xtc_filename (str) – Trajectory filename to read from and extend

• autocenter (bool) – Flag which, if set to True and there’s a single chain will center the protein in the box. This is useful for visualization purposes but does mean any translational diffusion will be lost. Default = False

Returns:

No return, but the existing XTC file is extended by one frame and then saved to disk.

Return type:

None

pimms.lattice_utils.build_all_envelope_pairs(positions, LR_binary_array, type_lattice, dimensions, hardwall=False)

Expects a LIST of positions and a numpy array of positions which engage in long-range interactions (or not) - 0 if not and 1 if yes.

Returns a list of tuples, where each tuple is a pair of positions.

The complete set of these positions represents the non-redudant set of long -range and short-range pairwise interactions associated with the positions defined in position

Parameters:

• positions (list) – list of positions

• LR_binary_array (numpy array) – numpy array of positions which engage in long-range interactions (or not) - 0 if not and 1 if yes.

• type_lattice (numpy.ndarray) – The type grid used by the inner-loop routines to determine which long-range / super-long-range pairs are generated from each position.

• dimensions (list) – A list of length 2 or 3 giving the lattice dimensions.

• hardwall (bool, optional) – If True, the hardwall inner-loop variants are used so that pairs do not straddle the periodic boundary. Default is False.

Returns:

A 3-tuple (SR_pairs, LR_pairs, SLR_pairs) of numpy arrays giving, respectively, the duplicate-free short-range, long-range and super-long-range interaction pairs. Each array has shape (N, 2, 2) in 2D or (N, 2, 3) in 3D. If positions is empty, three empty arrays of the appropriate shape are returned.

Return type:

tuple

pimms.lattice_utils.build_envelope_pairs(positions, dimensions, hardwall=False)

Expects a LIST of positions. Returns a unique unordered list of tuples, where each tuple is a pair of positions.

The complete set of these positions represents the non-redundant set of positions that make contact with the positions in the past $positions variable.

dimensions is the dimensions of the lattice

hardwall is a boolean which determines if we allow a pair of positions to straddle the boundary (in a periodic manner) or not

Parameters:

• positions (list) – A list of positions (each a 2D or 3D list) for which the enveloping short-range contact pairs are required.

• dimensions (list) – A list of length 2 or 3 giving the lattice dimensions.

• hardwall (bool, optional) – If True, pairs that straddle the periodic boundary are excluded (hardwall variant). Default is False.

Returns:

A numpy array of shape (N, 2, 2) in 2D or (N, 2, 3) in 3D, where each element is an unordered pair of positions making short-range contact with the input positions. Duplicate pairs are removed. An empty array of the appropriate shape is returned if positions is empty.

Return type:

numpy.ndarray

pimms.lattice_utils.center_of_mass_from_positions(positions, dimensions, on_lattice=True)

Return the center of mass from the list of positions. Assumes all positions have the same mass!

on_lattice can be set to True if you want a lattice-based COM or set to False if you want the true off-lattice Euclidean COM

COM is calculated by implementing the agorithm developed by Bai and Breen [1] extended to 3D, which means it determines the correct center of mass in a periodic box.

[1] Bai, L., & Breen, D. (2008). Calculating Center of Mass in an Unbounded 2D Environment. Journal of Graphics, GPU, and Game Tools, 13(4), 53 - 60.

Parameters:

• positions (list) – List of positions to calculate the center of mass from.

• dimensions (list) – List of dimensions of the box.

• on_lattice (bool) – If True, the center of mass is calculated on the lattice. If False, the center of mass is calculated in Euclidean space.

Returns:

The center of mass of the positions (2D or 3D list depending on if the input positions are 2D or 3D).

Return type:

list

pimms.lattice_utils.center_positions(positions, dimensions)

Returns the positions after centering them in the box. This is useful for visualisation and also for calculating the radial distribution function (RDF) as it ensures that the RDF is not biased by the positions being offset from the centre of the box.

Added in v0.1.34

Parameters:

• positions (list) – A list of lists, where each inner list is a position on the lattice.

• dimensions (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied, that reflects the lattice dimensions.

Returns:

A list of lists, where each inner list is a position on the lattice, that has been centered in the box.

Return type:

list

pimms.lattice_utils.check_all_chain_connectivity(list_of_chain_objects, dimensions, verbose=True)

Debugging function that takes a LATTICE.chains and LATTICE.dimensions pair from the simulation object to check the chain connectivity over all chains in the simulation.

Parameters:

• list_of_chain_objects (dict) – Dictionary of chain objects

• dimensions (list) – List of dimensions for the lattice

• verbose (bool) – Flag to print out debug information. Default = True

Returns:

No return, but will raise an error if any chain is not connected

Return type:

None

pimms.lattice_utils.check_chain_connectivity(chainID, chain_positions, dimensions, verbose=True)

Debugging function which ensures that a set of positions correspond to a valid, connected chain. Useful for debugging new moves, though not designed for performance during real simulations. If verbose is set to True, will print out a “CONNECTIVITY FINE” message assuming the function completes without issue.

Parameters:

• chainID (int) – Chain ID number

• chain_positions (list of lists) – List of lists of positions for the chain

• dimensions (list) – List of dimensions for the lattice

• verbose (bool) – Flag to print out debug information. Default = True

Returns:

No return, but will raise an error if the chain is not connected

Return type:

None

pimms.lattice_utils.close_xtc_writer(writer)

Close an open XTC writer (flushing the file). Safe to call with None.

Parameters:

writer (mdtraj.formats.XTCTrajectoryFile or None) – The writer handle to close.

Return type:

None

pimms.lattice_utils.convert_chain_to_single_image(chain_of_positions, dimensions)

If passed a chain of positions converts them so, relative to the first position, they’re all in the same image (single image convention). This assumes that each consecutive position is one lattice site apart from the next lattice site.

chain_of_positions should be a list of 2D or 3D positions, note that the first position in this list will define the PBC frame being used. We may want (in the future) to offer the option to use the terminal position, and then assess which reference frame is the least perturbative to the chain.

Also worth bearing in mind that if this method is applied to a set of chains in a connected cluster it will probably break everything. We have specific algorithm (snakesearch) to solve this problem, which can be found in the cluster_utils module.

dimensions should be list giving the dimensions of the box ([X,Y], or [X,Y,Z])

Parameters:

• chain_of_positions (list) – A list of 2D or 3D positions describing a chain, where each consecutive position is assumed to be one lattice site apart from the next. The first position defines the periodic image used as the reference frame.

• dimensions (list) – A list of length 2 or 3 giving the dimensions of the box, used as the box width in each dimension when unwrapping positions across periodic boundaries.

Returns:

A list of positions, the same length as chain_of_positions, unwrapped into a single periodic image and then offset so that all coordinates in every dimension are non-negative.

Return type:

list

Raises:

LatticeUtilsException – Raised if the chain cannot be unwrapped into a single image within the internal escape-counter limit, which indicates the input chain contains impossible (non-unit) bonds.

pimms.lattice_utils.delete_chain_by_ID(chainID, lattice_grid)

Deletes a chain from the lattice based on the chain’s ID

All lattice sites whose value equals chainID are reset to 0.0 (solvent).

Parameters:

• chainID (int) – The chain ID whose residues should be removed from the lattice.

• lattice_grid (numpy.ndarray) – The 2D or 3D lattice grid array, modified in place.

Return type:

None

pimms.lattice_utils.delete_chain_by_position(positions, lattice_grid, chainID=None)

Deletes a chain based on supplied position. Can be an entire chain or simply a portion of a chain

Parameters:

• positions (list) – A list of positions (each a 2D or 3D list) to be reset to solvent (0).

• lattice_grid (numpy.ndarray) – The 2D or 3D lattice grid array, modified in place.

• chainID (int, optional) – If provided, each position is checked to confirm it is currently occupied by this chain before being cleared; a mismatch raises an exception. If None (default) positions are cleared without checking.

Return type:

None

Raises:

ChainDeletionFailure – Raised (only when chainID is provided) if any position is not occupied by the expected chain.

pimms.lattice_utils.delete_residue(position, lattice, chainID=None)

Delete a residue at a given position in the lattice. This function will raise an exception if the position is already occupied by a residue from a different chain. This is the safe version of the function. If you want to overwrite the residue, set safe=False.

Parameters:

• position (list) – Position to delete the residue from.

• lattice (np.array) – Lattice to delete the residue from.

• chainID (int) – Chain ID of the residue to delete. If None, the residue will be deleted regardless of the chain ID.

Return type:

None

Raises:

ResidueAugmentException – Raised (only when chainID is provided) if the residue currently at position does not belong to the expected chain.

pimms.lattice_utils.do_positions_stradle_pbc_boundary(chain_positions)

For a set of positions returns true if the positions straddle a boundary else return False. Note this assumes the positions are connected to one another (i.e. each adjacent position is only 1 lattice site away from the next).

The check is performed by walking through consecutive positions and asking, in each dimension, whether the absolute difference between adjacent positions exceeds 1. A difference greater than 1 implies the bond wraps across a periodic boundary.

Parameters:

chain_positions (list) – A list of positions (each a length-2 or length-3 list) that are assumed to be consecutively connected along a chain.

Returns:

True if any consecutive pair of positions straddles a periodic boundary, otherwise False.

Return type:

bool

pimms.lattice_utils.find_nearest_position(target, positions_list, dimensions)

Given some target position (target) what is the index of the position in the positions list that is closest to that target? If multiple positions are found we simply return the first one in the positions list.

Parameters:

• target (list) – A 2D or 3D position (list of ints) to which all positions in positions_list are compared.

• positions_list (list) – A list of 2D or 3D positions (each a list of ints) to be compared against the target.

• dimensions (list) – A list of length 2 or 3 giving the lattice dimensions, used when computing the real (Euclidean, periodic) distance.

Returns:

A 2-tuple (int, float) where element [0] is the index of the position in positions_list closest to target, and element [1] is the actual distance between target and that closest position.

Return type:

tuple

Raises:

LatticeUtilsException – Raised if positions_list is empty.

pimms.lattice_utils.finish_pdb_file(filename)

Function that finalizes a PDB by adding terminating information.

Parameters:

filename (str) – Filename to be finalized

Returns:

No return but the file associated with filename is finalized as a PDB file.

Return type:

None

pimms.lattice_utils.get_adjacent_sites_2D(position1, position2, dimensions, extent_range=1)

Returns the lattice sites adjacent to a 2D position, generalized over an optional neighbourhood extent.

This is a thin wrapper around the compiled hyperloop.get_adjacent_sites_2D routine, which performs the periodic-boundary-corrected enumeration of neighbouring sites.

Parameters:

• position1 (int) – The x-coordinate of the position whose neighbours are required.

• position2 (int) – The y-coordinate of the position whose neighbours are required.

• dimensions (list) – A list of length 2 giving the lattice dimensions (X, Y).

• extent_range (int) – Half-width of the neighbourhood to enumerate around the position (default 1, i.e. immediately adjacent sites).

Returns:

A 3D numpy array where each element is a 2D numpy array describing an adjacent lattice position.

Return type:

numpy.ndarray

pimms.lattice_utils.get_adjacent_sites_3D(position1, position2, position3, dimensions, extent_range=1)

Returns the lattice sites adjacent to a 3D position, generalized over an optional neighbourhood extent.

This is a thin wrapper around the compiled hyperloop.get_adjacent_sites_3D routine, which performs the periodic-boundary-corrected enumeration of neighbouring sites.

Parameters:

• position1 (int) – The x-coordinate of the position whose neighbours are required.

• position2 (int) – The y-coordinate of the position whose neighbours are required.

• position3 (int) – The z-coordinate of the position whose neighbours are required.

• dimensions (list) – A list of length 3 giving the lattice dimensions (X, Y, Z).

• extent_range (int) – Half-width of the neighbourhood to enumerate around the position (default 1, i.e. immediately adjacent sites).

Returns:

A 4D numpy array where each element is a 3D numpy array describing an adjacent lattice position.

Return type:

numpy.ndarray

pimms.lattice_utils.get_all_chains_in_connected_component(chainID, lattice_grid, chainDict, threshold=None, useChains=True, hardwall=False)

Function which given a chainID, a dictionary of chain-to-position mappings, and a lattice grid will return the set of chains in the connected component containing chainID. Note a connected component is a heterotypic structure - i.e. we are looking for a connected component made up of any chains, not a single type of chain.

Parameters:

• chainID (int) – The chainID of the chain we initially are asking about

• lattice_grid (2D or 3D np.array) – Standard lattice grid

• chainDict (dictionary mapping chainIDs to a list of positions or to a chain object) – Dictionary containing a mapping of each chainID to either a list of positions associated with that chain, or the Chain object associated with that chainID

• threshold (int) – The max size of the connected component we are looking for. If None, this is ignored, but if set and we generate a connected component larger than this, we will raise a ClusterSizeThresholdException. Enables us to avoid situations where we’ve moving massive giant clusters around which may not be efficient if 90% of the chains are in the cluster.

useChainsBool

Boolean flag which defines if the chainDict is a true dictionary mapping chainID to a set of positions, or in fact a dictionary of Chain objects (which contain positions which must be accessed using the .get_ordered_positions()). This isn’t so much a feature as the fact that we want this function to be able to accept two different types of chain information (dictionary of lists of positions or dictionary of chain objects)

hardwallBool

Boolean flag which defines if we are using a hardwall potential or not. If we are, we will use a different method to calculate the connected component. If we are not, we will use a more efficient method which uses a union-find algorithm to calculate the connected

Returns:

A list of chainIDs associated with the chains in the connected component which contans the chain defined by $chainID

Return type:

list

pimms.lattice_utils.get_all_chains_in_long_range_cluster(chainID, latticeObject, hardwall=False)

Function which given a chainID, a dictionary of chain-to-position mappings, and a lattice grid will return the set of chains in the connected component where connectivity is defined in terms of long-range interactions. Note a connected component is a heterotypic structure - i.e. we are looking for a connected component made up of any chains, not a single type of chain.

Parameters:

• chainID (int) – The chainID of the chain we initially are asking about

• latticeObject (Lattice object) – Lattice object containing the lattice grid, the type grid, and the chain dictionary

• hardwall (Bool) – Boolean flag which defines if we are using a hardwall potential or not.

Returns:

A list of chainIDs associated with the chains in the connected component which contans the chain defined by $chainID

Return type:

list

pimms.lattice_utils.get_dimensions(lattice_grid)

Function that returns the dimensions associated with the lattice

Parameters:

lattice_grid (np.array) – A lattice grid numpy array

Returns:

Returns a numpy array of length 2 or 3, depending on the dimensionality of the system, where the value at each position reflects the size of the lattice in that dimension.

Return type:

np.array

pimms.lattice_utils.get_empty_site(lattice_grid, adjacentTo=None, hardwall=False)

Function which returns the position of an empty site on the lattice. If adjacentTo is populated then the returned site is adjacent to the the position defined by adjacentTo.

When adjacentTo is None the function performs random rejection sampling of the whole lattice until an empty site is found. When adjacentTo is set, only the sites neighbouring that position are considered (optionally excluding sites that straddle the boundary when hardwall is True), and one empty neighbour is selected at random.

Parameters:

• lattice_grid (numpy.ndarray) – The 2D or 3D lattice grid array, where a value of 0 denotes an empty (solvent) site.

• adjacentTo (list, optional) – If provided, a 2D or 3D position; the returned empty site is constrained to be adjacent to this position. If None (default), an empty site anywhere on the lattice is returned.

• hardwall (bool, optional) – If True (and adjacentTo is set), only neighbouring sites that do not straddle a periodic boundary are considered. Default is False.

Returns:

If adjacentTo is None, returns a single position (list) of an empty site. If adjacentTo is set, returns a 2-tuple (position, found) where position is the chosen empty neighbour (or a position filled with -1 if none was found) and found is a bool indicating success.

Return type:

list or tuple

Raises:

LatticeUtilsException – Raised (when adjacentTo is None) if the lattice is fully occupied or no empty site is found within the attempt limit.

pimms.lattice_utils.get_gridvalue(position, lattice_grid)

Returns the value on the lattice grid associated with the position defined by the 2/3 place tuple

The dimensionality (2D or 3D) is inferred from the shape of the lattice grid.

Parameters:

• position (list) – A 2D or 3D position (list of ints) to look up.

• lattice_grid (numpy.ndarray) – The 2D or 3D lattice grid array.

Returns:

The grid value stored at position (0 denotes an empty/solvent site, otherwise the occupying chain ID).

Return type:

int or float

Raises:

LatticeUtilsException – Raised if the lattice grid has an unsupported dimensionality.

pimms.lattice_utils.get_gridvalue_2D(position, lattice_grid)

Returns the value on the lattice grid associated with a 2D position.

This is a dimensionality-specialized variant of get_gridvalue() that assumes a 2D position and avoids the dimensionality check for speed.

Parameters:

• position (list) – A 2D position (list of two ints) to look up.

• lattice_grid (numpy.ndarray) – The 2D lattice grid array.

Returns:

The grid value stored at position (0 denotes an empty/solvent site, otherwise the occupying chain ID).

Return type:

int or float

pimms.lattice_utils.get_gridvalue_3D(position, lattice_grid)

Returns the value on the lattice grid associated with a 3D position.

This is a dimensionality-specialized variant of get_gridvalue() that delegates to the compiled hyperloop.get_gridvalue_3D routine for speed.

Parameters:

• position (list) – A 3D position (list of three ints) to look up.

• lattice_grid (numpy.ndarray) – The 3D lattice grid array.

Returns:

The grid value stored at position (0 denotes an empty/solvent site, otherwise the occupying chain ID).

Return type:

int or float

pimms.lattice_utils.get_real_distance(posA, posB, dimensions)

Function to calculate the real distance between two positions.

Parameters:

• posA (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied, that reflects a specific position on the lattice.

• posB (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied, that reflects a specific position on the lattice.

• dimensions (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied, that reflects the lattice dimensions.

Returns:

Returns a value that reflects the Euclidian distance between two positions on the lattice.

Return type:

float

pimms.lattice_utils.insert_chain(chainID, chain_length, lattice_grid, default_start=None, hardwall=False)

Function that inserts a chain into the passed lattice

Parameters:

• chainID (int) – Unique ID that identifies a specific chain on the lattice

• chain_length (int) – Number of residues in the chain

• lattice_grid (numpy.ndarray) – The 2D or 3D lattice grid array into which the chain is inserted. This array is modified in place as residues are placed.

• default_start (list, optional) – If provided, a position used as the fixed starting site for chain growth. If None (default) a random empty site is selected as the start.

• hardwall (bool, optional) – If True, chain growth only considers neighbouring sites that do not straddle a periodic boundary. Default is False.

Returns:

A list of positions describing the inserted chain, ordered along the chain. The lattice grid is also updated in place.

Return type:

list

Raises:

ChainInsertionFailure – Raised if a valid chain configuration could not be constructed within CONFIG.CHAIN_INIT_ATTEMPTS attempts.

pimms.lattice_utils.insert_residue(position, lattice, chainID, safe=True)

Insert a residue at a given position in the lattice. This function will raise an exception if the position is already occupied by a residue from a different chain. This is the safe version of the function. If you want to overwrite the residue, set safe=False.

Parameters:

• position (list) – Position to insert the residue

• lattice (np.array) – Lattice to insert the residue into

• chainID (int) – Chain ID to insert the residue for

• safe (bool) – If True, will raise an exception if the position is already occupied. If False, will overwrite the residue.

Return type:

None

pimms.lattice_utils.make_chain_whole(chain_of_positions, dimensions)

Unwrap a chain into a single periodic image, anchored at the FIRST bead’s real (in-box) position.

Each consecutive bead is placed exactly one lattice step from the previous one, so the chain is never torn across a periodic boundary; coordinates may fall outside the box on either side. Unlike convert_chain_to_single_image(), this does NOT afterwards shift the whole chain to be non-negative - that shift would translate every boundary-crossing chain toward one face of the box, which is why an unwrapped trajectory built with the single-image routine only ever appeared to bulge out of one side. Keeping the first bead in place makes chains spill symmetrically out of whichever face they actually cross.

This is used purely for trajectory visualisation (TRAJECTORY_PBC_UNWRAP) and is a lightweight, allocation-cheap loop (no deep copies / numpy transposes), since it is called once per chain per written frame.

Parameters:

• chain_of_positions (list) – List of 2D or 3D integer positions describing a chain, consecutive beads one lattice site apart. The first position defines the reference image.

• dimensions (list) – Box dimensions (length 2 or 3), used as the per-axis period.

Returns:

The chain positions unwrapped into a single image, first bead unchanged.

Return type:

list

Raises:

LatticeUtilsException – If a bond cannot be resolved within the escape-counter limit (indicates an impossible / non-unit bond in the input chain).

pimms.lattice_utils.open_pdb_file(dimensions, spacing, filename='lattice.pdb')

Function that initializes a PDB file to be written to.

Parameters:

• dimensions (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied, that reflects the lattice dimensions.

• spacing (float) – Lattice-to-realspace spacing in angstroms.

• filename (str) – Filename to write to

Returns:

No return value, but a new PDB file is initialized on disk.

Return type:

None

pimms.lattice_utils.open_xtc_writer(lattice, spacing, pdb_filename='START.pdb', xtc_filename='traj.xtc', autocenter=False, unwrap=False)

Write the topology PDB, open a persistent XTC writer, write the first frame, and return the open writer handle.

This is the efficient replacement for the previous per-frame append_to_xtc_file_non_redundant approach, which re-loaded the entire growing trajectory from disk and re-saved it on EVERY frame - O(frames^2) in both wall time and disk I/O. Here a single mdtraj XTC file handle is kept open for the whole run and each frame is appended with write_xtc_frame() in O(1); the handle is closed with close_xtc_writer().

Parameters:

• lattice (lattice.Lattice) – Current lattice.

• spacing (float) – Lattice-to-realspace spacing (angstroms).

• pdb_filename (str) – Topology PDB filename to (re)write.

• xtc_filename (str) – Trajectory filename to create.

• autocenter (bool) – Single-chain autocentring (see build_pdb_file). Default False.

• unwrap (bool) – Make chains whole across PBC before writing (TRAJECTORY_PBC_UNWRAP). Default False.

Returns:

The open writer handle (write more frames with write_xtc_frame, then close with close_xtc_writer).

Return type:

mdtraj.formats.XTCTrajectoryFile

pimms.lattice_utils.pbc_convert(position, dimensions)

Returns lattice site positions after carrying out periodic boundary conditions (PBC) conversions.

Parameters:

• position (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied, that reflects a specific position on the lattice.

• dimensions (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied, that reflects the lattice dimensions used as the modulus for each dimension.

Returns:

A list where the length matches the input position list, where the new positions reflects the periodic-boundary condition corrected positions.

Return type:

list

pimms.lattice_utils.pbc_correct(posA, posB, dimensions)

Returns the two positions after converting them such that, relative to one another, they are in the same image (i.e. not cross a border). NOTE that this function ONLY allows works for a single periodic image but if you had some weird set of positions that are spanning multiple images this won’t work. This is a pretty unlikely but just FYI.

For simplicity the method assumes posA is God and posB can be re- set. This is abirtrary, but we don’t need to futz with them both!!

NOTE: This method was introduced in 0.9 so has not been as throughly vetted as a lot of the core code in PIMMS.

Parameters:

• posA (list) – A list of length 2 or 3, depending on the dimensionality of the system, giving a position on the lattice. This position is treated as fixed.

• posB (list) – A list of length 2 or 3 (matching posA) giving a second position on the lattice. This position may be shifted by a single box length in any dimension so that it lies in the same periodic image as posA.

• dimensions (list) – A list of length 2 or 3 giving the lattice dimensions, used as the box width in each dimension.

Returns:

A 2-tuple (posA, newB) where posA is the unchanged input position and newB is posB shifted (where required) so that, in each dimension, the separation between the two positions is the minimum image separation.

Return type:

tuple

pimms.lattice_utils.place_chain_by_position(positions, lattice_grid, chainID, safe=False)

Sets the positions defined in the positions vectors to a chain. Note this can be an entire chain or a subset of a chain

Parameters:

• positions (list) – A list of positions (each a 2D or 3D list) that will be assigned to the chain.

• lattice_grid (numpy.ndarray) – The 2D or 3D lattice grid array, modified in place.

• chainID (int) – The chain ID value written into the lattice grid at each position.

• safe (bool, optional) – If True, each target position is checked to ensure it is currently empty before writing; an occupied site raises an exception. If False (default) positions are overwritten without checking.

Return type:

None

Raises:

ChainInsertionFailure – Raised (only when safe is True) if any target position is already occupied by another chain.

pimms.lattice_utils.rotate_positions_2D(positions, degrees)

Functions to carry out 2D cardinal position rotation around the origin.

The CARDINAL_ROTATION_2D matrix is assigned in CONFIG, affording extremely fast rotation.

Parameters:

• positions (list) – A list of 2D positions to be rotated.

• degrees (int) – The number of degrees to rotate the positions by. Must be one of 90, 180, or 270.

Returns:

rotated_positions – A list of 2D positions that have been rotated by the specified number of degrees.

Return type:

list

pimms.lattice_utils.rotate_positions_3D(positions, dimension, degrees)

Functions to carry out cardinal position rotation around the origin.

The CARDINAL_ROTATION_3D matrix is assigned in CONFIG, affording extremely fast rotation.

Parameters:

• positions (list) – List of positions to rotate

• dimension (str) – Dimension to rotate around. Must be one of ‘x’, ‘y’, or ‘z’.

• degrees (int) – Degrees to rotate by. Must be one of 90, 180, or 270.

Returns:

List of rotated positions

Return type:

list

pimms.lattice_utils.run_rotation(positions, rotation_matrix)

low-level function that performs single point rotation. This should generally not be called but instead the wrapper functions rotate_positions_3D or rotate_positions_2D should be used.

Parameters:

• positions (list) – List of positions to rotate

• rotation_matrix (np.array) – Rotation matrix to apply

Returns:

List of rotated positions

Return type:

list

pimms.lattice_utils.same_sites(site1, site2)

Explicit function to check two sites are the same. Automatically generalizes to 2D or 3D in response to the site dimensions.

Parameters:

• site1 (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied.

• site2 (list) – A list of length 2 or 3, depending on the dimensionality of the system being studied.

Returns:

Returns True if the two sites are the same and False if not

Return type:

bool

pimms.lattice_utils.save_out_sim(master_traj, xtc_filename)

Save out the master trajectory.

Parameters:

• master_traj (mdtraj.Trajectory) – master trajectory we will build throught the sim.

• xtc_filename (str) – Filename to write to disk.

Returns:

No return, but the existing XTC file is saved to disk

Return type:

None

pimms.lattice_utils.set_gridvalue(position, value, lattice_grid)

Sets the position defined at lattice site $position on $lattice_grid to $value. This CHANGES the $lattice_grid object (which is assumed to be a numpy 2D or 3D array) and returns it.

Parameters:

• position (list) – A 2D or 3D position (list of ints) at which to write.

• value (int or float) – The value to write into the lattice grid at position (e.g. a chain ID or 0 for solvent).

• lattice_grid (numpy.ndarray) – The 2D or 3D lattice grid array, modified in place.

Returns:

The same lattice grid object that was passed in, after modification.

Return type:

numpy.ndarray

Raises:

LatticeUtilsException – Raised if position has an unsupported dimensionality.

pimms.lattice_utils.start_xtc_file(lattice, spacing, pdb_filename='START.pdb', xtc_filename='traj.xtc', unwrap=False)

Function that initializes a new .xtc file. This deletes an existing XTC file of the same name to avoid any issues.

Parameters:

• lattice (lattice.Lattice) – Current Lattice object

• spacing (float) – Lattice-to-realspace spacing in angstroms, used when writing the corresponding PDB file.

• pdb_filename (str) – New XTC files need a corresponding PDB file. This defines the name of that PDB file.

• xtc_filename (str) – New XTC files need a corresponding PDB file. This defines the name of that PDB file.

Returns:

No return value, but a newly initialized XTC file is generated

Return type:

None

pimms.lattice_utils.update_master_traj(lattice, spacing, master_traj, pdb_filename, autocenter=False, unwrap=False)

Low level function that adds a current lattice to an existing XTC file. This is different than ‘append_to_xtc_file’ in that it does not read in an existing XTC file but instead appends a frame to the passed master trajectory object.

If the master_traj object has not yet been initialized, this will also read in the pdb_filename and initialize the master_traj object using that as a topology file. This means PDB initialization needs to happen BEFORE this function is called! This is a deliberate design choice to allow for the master_traj object to be passed between functions and have the topology file be read in only once.

Parameters:

• lattice (lattice.Lattice) – A lattice object

• spacing (float) – Lattice-to-realspace spacing in angstroms.

• master_traj (mdtraj.Trajectory) – master trajectory we will build throught the sim.

• pdb_file_name (current_pdb_filename) – the current_pdb_filename from simulation.py

• autocenter (bool) – Flag which, if set to True and there’s a single chain will center the protein in the box. This is useful for visualization purposes but does mean any translational diffusion will be lost. Default = False

Returns:

Returns the master trajectory object after update, although because this master trajectory is passed by value in principle this return object does not need to be dealt with as the passed object is updated in place.

Return type:

mdtraj.Trajectory

pimms.lattice_utils.write_lattice_to_pdb(latticeObject, spacing, filename='lattice.pdb', write_connect=False, autocenter=False, unwrap=False)

Wrapper function that dumps the current Lattice object to a PDB file

Parameters:

• latticeObject (lattice.Lattice) – Current lattice object

• spacing (float) – Lattice-to-realspace spacing in angstroms.

• filename (str) – Filename to write to (default is lattice.pdb)

• write_connect (bool) – Flag to write connect information to the PDB file (default is False)

• autocenter (bool) – Flag to center the lattice in the PDB file (default is False). Note that this correctly is dealth with in build_pdb_file - if more than one chain this is ignored.

Return type:

None

pimms.lattice_utils.write_xtc_frame(writer, lattice, spacing, autocenter=False, unwrap=False)

Append one frame from the current lattice to an open XTC writer (O(1), no reload).

Parameters:

• writer (mdtraj.formats.XTCTrajectoryFile) – Open writer handle from open_xtc_writer().

• lattice (lattice.Lattice) – Current lattice.

• spacing (float) – Lattice-to-realspace spacing (angstroms).

• autocenter (bool) – Single-chain autocentring. Default False.

• unwrap (bool) – Make chains whole across PBC before writing. Default False.

Return type:

None

Input parsing & configuration

class pimms.keyfile_parser.KeyFileParser(filename, parse_only=False)

Bases: object

KeyFileParser is essentially where all the logic that deals with input information is defined.

Specifically, a KeyFileParser object can read a keyfile and extract any/all pertinant information there in. It then sets any default values which can be estimated but weren’t defined explicitly. Finally, it sanity checks all the input to ensure we’re doing something sensible.

add_derived_keywords()

Add keywords that are derived from already-parsed keyword values.

This runs after parsing, default assignment and sanity checking. It sets the end-to-end-distance analysis frequency (ANA_END_TO_END) equal to the polymer-analysis frequency, and sets EQUILIBRIUM_TEMPERATURE to the final simulation temperature (the quench end temperature for a quench run, otherwise the standard TEMPERATURE).

Returns:

No return value, but self.keyword_lookup is updated in place with the derived keywords ANA_END_TO_END and EQUILIBRIUM_TEMPERATURE.

Return type:

None

assign_default()

Build the self.DEFAULTS dictionary of default keyword values.

The absolute reference default values are taken from CONFIG.DEFAULTS. In reality these could live in a separate configuration file, but in the interest of keeping everything within this file we define those default values here.

This basically builds up a self.DEFAULTS dictionary which defines default values for each keyword. Note that for a few of these the default can depend on the keywords parsed.

Note that these defaults are used in two ways:

1. If keywords are not provided then they are used as the default. Obviously.

2. This can be used as a test to ask if a keyword has been set, because IF you pass a keyword and it changes the parsed keyword from the default this is the only time we care about a keyword being provided, hence functionaly this is how we define if a keyword is provided (or not). In particular, this is used for evaluating if EXPERIMENTAL keywords are being used.

Returns:

No return value, but self.DEFAULTS is populated in place (including analysis-frequency defaults derived from ANALYSIS_FREQ and a freshly generated random SEED).

Return type:

None

parse(filename, verbose=True)

Main function reads in a keyfile and extracts out the relevant details based on the keywords. Keywoerds are assigned to the self.keywords_lookup

Keywords must be defined as

KEYWORD : VALUE # comment

This function reads through all the lines in the keyfile and in a keyword-specific manner parses each keyword and assigns it to the self.keyword_lookup dictionary. Any keywords that are expected but NOT defined in the keyfile are then assigned as default values later.

Parameters:

• filename (str) – Name of the keyfile to be read

• verbose (bool (default = True)) – Flag which determines the level of warning messages to print.

Returns:

No return type but assigns values to the self.keyword_lookup dictionary, which itself is a key-value pair fo simulations keywords

Return type:

None

print_summary()

Print a full human-readable summary of the parsed simulation setup.

Reports the system overview (step counts, temperatures, expected number of frames), box dimensions and resulting occupied volume fraction / solute concentration, quench settings, freeze-file settings, simulation components (chains), and output frequencies. A subset of these values are also written to the PIMMS logfile via pimmslogger.

Returns:

No return value; the summary is printed to standard output (and key values are logged).

Return type:

None

run_sanity_checks()

Function which, once the keyword has been parsed, will run an arbitrary number of sanity checks to ensure things make sense. Adding new checks is simply a case of adding a new ‘CHECK’ section in the code.

This also sets a number of ‘internal’ keywords - keywords which are derived from the manually defined keywords but are NOT explicitly read in. This lets the system set certain configurations by deducing things from the manual keyfile that aid in decision making. To keep things clear a list of the ‘internal’ keywords are below. Internal keywords always begin with a ‘__’

__TSMMC_USED # BOOLEAN –> set to true or false depending on if any TSMMC are being used.

Specific details of each sanity check are provided as code blocks below. Before adding additional checks please read through these! Note also that restart file sanity checks are done in their own function which gets called from in here

Returns:

No return value, but self.keyword_lookup is updated in place (e.g. derived __ keywords are set, chains may be upper-cased, restart/freeze files are loaded, and quench parameters normalised).

Return type:

None

Raises:

KeyFileException – If any sanity check fails (e.g. out-of-range numerical values, move fractions that do not sum to 1.0, missing/invalid parameter or restart files, or incompatible box/hardwall/experimental settings).

sanity_check_and_update_with_restart_file()

Run sanity checks and update the CHAIN information so the box dimensions/concentration info is correct. Note that we DON’T sanity check the restart file input against the keyfile-defined chain values. This means we can be agnostic about what’s in a restart file, and allows restarts to be more permissive.

This set of sanity checks is actually the last thing done, which means that we have to re-check some of the things that were already checked with updated information read from the restart file.

Below is the general ruberic for how restart files are dealt with:

• The restart file fully overwrites all chain information. Chain information in a keyfile is completely ignored if a restart file is provided.

• By default we assume the keyfile provided dimensions and hardwall status are to be used by in the restart. However, this may actually not be compatible, in which case an exception is thrown. The RESTART_OVERRIDE_DIMENSIONS, RESTART_OVERRIDE_HARDWALL force the simulation to use the dimensions and hardwall values passed by the restart file. The one trick here is that if the keyfile requires an resized_equilibration then the restart file’s dimensions will be used for the initial equilibration, and the production part of the simulation will be run using information from the keyfile.

Returns:

No return value. self.keyword_lookup is updated in place (the CHAIN composition is rebuilt from the restart object, and HARDWALL / DIMENSIONS may be overridden), and the restart object’s lattice dimensions may be recentred/resized. Returns early (without changes) if no restart file is set.

Return type:

None

Raises:

RestartException – If the restart file is incompatible with the keyfile (e.g. mismatched dimensionality, conflicting hardwall/PBC modes, or box dimensions that are smaller than the restart file’s dimensions).

set_defaults()

Assign default values from self.DEFAULTS to any unset keywords.

First the presence of every required keyword is validated (and at least one of CHAIN or RESTART_FILE must be present). Then every expected keyword that was not explicitly supplied in the keyfile is set to its default value from self.DEFAULTS, announcing each default as it is applied.

Returns:

No return value, but self.keyword_lookup is updated in place with default values for any keywords that were not explicitly defined.

Return type:

None

Raises:

KeyFileException – If a required keyword is missing, or if neither CHAIN nor RESTART_FILE was provided.

set_dynamic_defaults()

This function allows final default values to update such that ‘defaults’ can respond to passed parameters.

Parameters:

No parameters, but updates the self.keyword_lookup dictionary with new default values as needed

Return type:

No return value, but updates the self.keyword_lookup dictionary with new default values as needed

update_keyfile(update_dictionary)

Function that takes an update dictionary of key : value pairs and updates the current key:value pairs with these new keywords.

Parameters:

update_dictionary (dictionary) – key-value where keys are keywords and values are the values to change to. Any key value pair will be added

Return type:

No return value, but updates the self.keyword_lookup dictionary

write_keyfile(output_filename, PADDING=10)

Write the key-value pairs from the keyword_lookup dictionary to a file.

Each keyword and its value are written on a single line in KEY : value format, with padding inserted between the key and the value for readability.

Parameters:

• output_filename (str) – The name (path) of the file to write the key-value pairs to.

• PADDING (int, optional) – The minimum number of spaces to use for padding between the key and value. Default is 10.

Returns:

No return value; the keyword/value pairs are written to disk.

Return type:

None

Examples

>>> parser = KeyFileParser('input.kf')
>>> parser.write_keyfile('output.txt')

pimms.keyfile_parser.print_keyword_info()

Print a human-readable table of every supported keyfile keyword.

Iterates over CONFIG.KEYWORDS_DESCRIPTION and prints, for each keyword, its name, its expected value type, and a short description. Output is written to standard output in aligned columns.

Parameters:

None

Returns:

No return value; the keyword table is printed to standard output.

Return type:

None

pimms.parameterfile_parser.parse_angles(filename, temperature=False)

Reads in a parameter file and constructs an angle penalty dictionary. We allow the definition of two types of angle penalties

ANGLE_PENALTY lines define ABSOLUTE (integer) angle-penalty values, while ANGLE_PENALTY_T_NORM lines define temperature-normalised (kT, assuming k = 1) penalties that are multiplied by temperature at parse time to give the absolute penalty. All non-angle lines in the file are ignored.

Parameters:

• filename (str) – Filename for the parameter file to read.

• temperature (float or bool, optional) – Simulation temperature used to scale ANGLE_PENALTY_T_NORM entries. Must be a numeric value if any T-normalised lines are present; defaults to False (i.e. unset).

Returns:

Dictionary keyed by residue name, mapping each residue to a list of three angle-penalty values [AP1, AP2, AP3].

Return type:

dict

Raises:

ParameterFileException – If an angle-penalty line is malformatted, defines non-numeric penalty values, redefines a residue that already has a penalty, or uses ANGLE_PENALTY_T_NORM without a temperature being supplied.

pimms.parameterfile_parser.parse_energy(filename)

Function which reads in an energy parameter file and returns the interaction matrix as a redundant dictionary of dictionaries and the set of non-redundant particles defined therein.

Each line in the parameter file defines the interaction between two residues, in the format

A B X

or

A B X Y

In both cases X defines the short range interaction energy between A and B. This occurs when A and B are adjacent to one another on a lattice

Y defines the long range interaction energy (i.e. electrostatics though you could use it for anything) which occurs over one lattice site.

Parameters:

filename (str) – Filename for the parameter file.

Returns:

A 5-tuple (energy_pairs, residue_names, LR_energy_pairs, LR_residue_names, SLR_energy_pairs) where:

• energy_pairsdict of dict

  Fully redundant short-range interaction matrix keyed energy_pairs[P1][P2] -> int.

• residue_nameslist of str

  Non-redundant residue type names, with solvent ('0') guaranteed to be the first element.

• LR_energy_pairsdict of dict

  Fully redundant long-range interaction matrix (missing pairs filled with 0.0).

• LR_residue_nameslist of str

  Residue type names that participate in long-range interactions.

• SLR_energy_pairsdict of dict

  Fully redundant semi-long-range interaction matrix, with the same key structure as LR_energy_pairs.

As a side effect, a copy of the parsed parameter file is written to CONFIG.OUTPUT_USED_PARAMETER_FILE for reproducibility.

Return type:

tuple

Raises:

ParameterFileException – If a line is malformed, a float is used as an interaction strength, the matrix is redundant or incomplete, no solvation interactions are defined, a non-zero solvent-solvent interaction is given, or long-range solvent-solute interactions are defined.

pimms.parameterfile_parser.write_angle_parameter_summary(angle_dict, filename)

Write a human-readable summary of the angle penalties to disk.

The summary is written to CONFIG.OUTPUT_FULL_ANGLE_POTENTIAL and lists, for each residue, its three angle-penalty values.

Parameters:

• angle_dict (dict) – Dictionary keyed by residue name mapping to a list of three angle-penalty values [AP1, AP2, AP3] (as produced by parse_angles()).

• filename (str) – Name of the source parameter file. Currently accepted for interface consistency but not written into the output.

Returns:

No return value; the angle-penalty summary is written to disk.

Return type:

None

class pimms.restart.RestartObject

Bases: object

Object used to read and write restart files. Restart information ONLY contains information on chain position, sequence, and type, and grid dimenisons, but does NOT include any information

Note that the self.chains object in a RestartObject has the following structure:

1. Is a dictionary

2. Keys are chainID (i.e. each seperate chain has it’s own entry)

3. values is a list with three elements [0] : bead positions (N->C) [1] : chain sequence (which will be referenced against the parameter file) [2] : chainType : a single value that defines the type of chain

add_extra_chains(extra_chains, log=False)

Function which allows extra chains (as read from a keyfile) to be added to a RestartObject so that when a new lattice is initialized from this RestartObject those extra chains are randomly placed somewhere across the simulation box.

Note extra_chains ONLY have a sequence and chainType associated with them, but do NOT have any positions.

Parameters:

• extra_chains (list) – List with two elements [0] = number of chains (int) [1] = chain sequence (str)

• log (bool) – Flag which, if set to true, means if this seq already has a chainType defined but the passed chainType is a DIFFERENT value it’ll warn the user about this.

Returns:

No return type, but the internal self.extra_chains dictionary will be appropriately updated.

Return type:

None

build_from_file(filename, log=False)

Function that constructs a restart object from a passed filename. Performs some sanity check in reading in the file but doesn’t actually check that the chain positions make sense on the lattice. We can and should probably make this better going forwards…

Parameters:

• filename (str) – Name of the file to be read

• log (bool (default = False)) – Flag which if set to True means warnings are written to the standard PIMMS logfile

Returns:

• None but updates the current object to contain self.dimensions, self.energy, self.hardwall

• and self.chains[] info.

build_from_lattice(LATTICE, hardwall=False, log=False)

Construct a restart object using a lattice object to set the chain positions.

Parameter

LATTICEpimms.lattice.Lattice

A standard PIMMS latticd object

hardwallbool (default = False)

Flag which sets of the current system defines a hardwall or, if false PBC.

logbool (default = False)

Flag which if set to True means warnings are written to the standard PIMMS logfile

returns:

No return type, but the internal self.chains dictionary will be appropriately updated.

rtype:

None

set_energy(energy)

Set the RestartObject’s stored energy value.

Parameters:

energy (float) – The system energy to record in the restart object (written out when the restart file is saved).

Returns:

No return value, but self.energy is updated in place.

Return type:

None

update_lattice_dimensions(new_dimensions, manual_offset=None)

Resize the lattice and reposition the chains within the new lattice.

Updates the restart object’s dimensions and shifts every chain position by a per-dimension offset. By default the offset is computed so the existing chains end up centred in the larger lattice; alternatively an explicit manual_offset can be supplied. If the offset would move any bead outside the new lattice, the original dimensions are restored and the underlying RestartException is re-raised.

Parameters:

• new_dimensions (list of int) – The new lattice dimensions (length 2 or 3). Should be greater than or equal to the current dimensions for centring to make sense.

• manual_offset (list of int, optional) – Explicit per-dimension offset to apply to every chain position. If None (the default), a centring offset is computed automatically from the difference between new_dimensions and the current dimensions.

Returns:

No return value, but self.dimensions and the stored chain positions are updated in place.

Return type:

None

Raises:

RestartException – If applying the offset would place a bead outside the new lattice (in which case the prior dimensions are restored before re-raising).

write_to_file()

Serialize the restart object to disk as a pickle file.

Writes a dictionary containing the chain information (CHAINS), lattice dimensions (DIMENSIONS), recorded energy (ENERGY) and hardwall flag (HARDWALL) to CONFIG.RESTART_FILENAME using pickle. Note that extra_chains are not written; only the materialised self.chains are saved.

Returns:

No return value; the restart data is written to CONFIG.RESTART_FILENAME.

Return type:

None

class pimms.data_structures.AnalysisSettings(cluster_threshold)

Bases: object

Lightweight container for on-the-fly analysis configuration.

Holds settings used by the analysis machinery during a simulation, currently just the cluster-size threshold.

class pimms.data_structures.FreezeFile(filename)

Bases: object

Reader/container for a simulation freeze file.

Parses a freeze file (see __init__() for the file format) and exposes the set of frozen chain IDs (and, in future, bead IDs) so that the simulation can hold those chains fixed during sampling.

property beads

The deduplicated bead IDs to be frozen (bead-level freezing is not yet implemented, so this is typically empty).

Type:

list of int

property chains

The deduplicated chain IDs to be frozen.

Type:

list of int

property filename

The path of the freeze file that was read.

Type:

str

log_freeze_file()

Write a summary of the freeze file to the simulation log.

Logs the freeze file path, the number of frozen chains, and the list of frozen chain IDs as STATUS entries.

Returns:

No return value; entries are appended to the simulation log.

Return type:

None

validate_freeze_file(latticeObject)

Function to validate that the chains and beads specified in the freeze file are actually present in the lattice object.

Parameters:

latticeObject (Lattice) – The lattice object to be validated against

Returns:

No return variable, but an exception is raised if the freeze file is not valid

Return type:

None

Analysis & output

pimms.analysis_IO.write_LR_cluster_properties(step, LR_cluster_polymeric_properties_list, LR_cluster_size_list, LR_cluster_radial_density)

Write the per-cluster polymeric and gross properties for long-range clusters.

Identical in behaviour to write_cluster_properties() but writes to the long-range (LR) cluster property output files.

Parameters:

• step (int) – Current simulation step.

• LR_cluster_polymeric_properties_list (list of list) – One entry per LR cluster; each entry is [Rg, asphericity].

• LR_cluster_size_list (list of list) – One entry per LR cluster; each entry is [volume, surface_area, density].

• LR_cluster_radial_density (list of list of float) – One entry per LR cluster; each entry is the radial density profile.

Returns:

Nothing is returned; results are appended to the LR cluster property files defined in CONFIG (OUTNAME_LR_CLUSTER_RG, OUTNAME_LR_CLUSTER_ASPH, OUTNAME_LR_CLUSTER_VOL, OUTNAME_LR_CLUSTER_AREA, OUTNAME_LR_CLUSTER_DENSITY and OUTNAME_LR_CLUSTER_RADIAL_DENSITY_PROFILE).

Return type:

None

pimms.analysis_IO.write_LR_clusters(step, clusters, IDtoType)

Write long-range (LR) cluster size distribution and composition for a step.

Identical in behaviour to write_clusters() but writes to the long-range cluster output files. LR clusters are defined as clusters connected through short-range OR long-range interactions.

Parameters:

• step (int) – Current simulation step.

• clusters (list of list of int) – List where each sublist represents a cluster containing the chainIDs of its member chains. For example [[1, 2], [3], [4]] represents a four-chain system where chains 1 and 2 are clustered together while chains 3 and 4 are isolated.

• IDtoType (dict) – Dictionary mapping each chainID to its chainType.

Returns:

Nothing is returned; results are appended to the LR cluster output files defined in CONFIG (OUTNAME_LR_CLUSTERS, OUTNAME_NUM_LR_CLUSTERS and per-chain-type CHAIN_<type>_ prefixed files).

Return type:

None

pimms.analysis_IO.write_acceptance_statistics(step, acceptanceObject)

Write move attempt, acceptance, and total move statistics for a step.

Writes the per-move attempt counts and accepted counts to two separate files (so attempted and accepted moves can be compared directly), and the cumulative total number of attempted MC moves (across all sub-loops) to a third file.

Parameters:

• step (int) – Current simulation step.

• acceptanceObject (AcceptanceCalculator) – Object tracking move statistics. Must expose move_count and accepted_count (indexable by move code) and alt_Markov_chain_moves.

Returns:

Nothing is returned; lines are appended to the files defined by CONFIG.OUTNAME_MOVES, CONFIG.OUTNAME_ACCEPTANCE and CONFIG.OUTNAME_TOTAL_MOVES.

Return type:

None

Raises:

AcceptanceException – If the number of tracked moves is not the expected hard-coded value (15), which guards the explicit move list used to compute total moves.

pimms.analysis_IO.write_asphericity(step, asphericity_list)

Append the instantaneous asphericity values for a single step.

Parameters:

• step (int) – Current simulation step.

• asphericity_list (list of float) – Asphericity for each chain (or analysed object) at this step.

Returns:

Nothing is returned; a line is appended to the file defined by CONFIG.OUTNAME_ASPH.

Return type:

None

pimms.analysis_IO.write_cluster_properties(step, cluster_polymeric_properties_list, cluster_size_list, cluster_radial_density)

Write the per-cluster polymeric and gross properties for a single step.

Appends one line per output file for the current step, distributing the supplied per-cluster properties across the dedicated short-range cluster property files (Rg, asphericity, volume, surface area, density and radial density profile).

Parameters:

• step (int) – Current simulation step.

• cluster_polymeric_properties_list (list of list) – One entry per cluster; each entry is [Rg, asphericity] as returned by the polymeric-property analysis.

• cluster_size_list (list of list) – One entry per cluster; each entry is [volume, surface_area, density] as returned by the gross-property analysis.

• cluster_radial_density (list of list of float) – One entry per cluster; each entry is the radial density profile (a list of densities as a function of distance from the cluster COM).

Returns:

Nothing is returned; results are appended to the cluster property files defined in CONFIG (OUTNAME_CLUSTER_RG, OUTNAME_CLUSTER_ASPH, OUTNAME_CLUSTER_VOL, OUTNAME_CLUSTER_AREA, OUTNAME_CLUSTER_DENSITY and OUTNAME_CLUSTER_RADIAL_DENSITY_PROFILE).

Return type:

None

pimms.analysis_IO.write_clusters(step, clusters, IDtoType)

Write short-range cluster size distribution and composition for a step.

Appends (1) the number of chains in each cluster, (2) the total number of clusters, and - when more than one chain type is present - (3) one file per chain type giving the fraction of each cluster made up of that chain type.

Parameters:

• step (int) – Current simulation step.

• clusters (list of list of int) – List where each sublist represents a cluster containing the chainIDs of its member chains. For example [[1, 2], [3], [4]] represents a four-chain system where chains 1 and 2 are clustered together while chains 3 and 4 are isolated.

• IDtoType (dict) – Dictionary mapping each chainID to its chainType.

Returns:

Nothing is returned; results are appended to the cluster output files defined in CONFIG (OUTNAME_CLUSTERS, OUTNAME_NUM_CLUSTERS and per-chain-type CHAIN_<type>_ prefixed files).

Return type:

None

pimms.analysis_IO.write_distance_map(dMap, prefix=False)

Write a square inter-residue distance map to file.

Writes (overwriting any existing file) the full seqlen x seqlen matrix as tab-separated rows, one row per source residue.

Parameters:

• dMap (numpy.ndarray) – Square (seqlen, seqlen) distance map matrix.

• prefix (str or bool, optional) – If False (default) the standard output filename is used. If a string is supplied it is prepended to the output basename.

Returns:

Nothing is returned; results are written to the file defined by CONFIG.OUTNAME_DMAP (optionally prefixed).

Return type:

None

pimms.analysis_IO.write_end_to_end(step, e2e_list)

Append the instantaneous end-to-end distance values for a single step.

Parameters:

• step (int) – Current simulation step.

• e2e_list (list of float) – End-to-end distance for each chain at this step.

Returns:

Nothing is returned; a line is appended to the file defined by CONFIG.OUTNAME_E2E.

Return type:

None

pimms.analysis_IO.write_energy(step, energy)

Append the current step and system energy to the ENERGY output file.

Parameters:

• step (int) – Current simulation step.

• energy (float) – Current total system energy.

Returns:

Nothing is returned; a line is appended to the file defined by CONFIG.OUTNAME_ENERGY.

Return type:

None

pimms.analysis_IO.write_internal_scaling(mean_IS, mean_IS_squared, prefix=False)

Write the mean internal scaling and mean internal scaling squared profiles.

Each profile is written (overwriting any existing file) as a two-column table of gap (sequence separation, starting at 1) versus the mean value at that gap.

Parameters:

• mean_IS (list of float) – Mean internal scaling value for each sequence-separation gap, ordered by increasing gap.

• mean_IS_squared (list of float) – Mean internal scaling squared value for each sequence-separation gap, ordered by increasing gap.

• prefix (str or bool, optional) – If False (default) the standard output filenames are used. If a string is supplied it is prepended to the output basenames (e.g. to separate per-chain-type output).

Returns:

Nothing is returned; results are written to the files defined by CONFIG.OUTNAME_INTERNAL_SCALING and CONFIG.OUTNAME_INTERNAL_SCALING_SQUARED (optionally prefixed).

Return type:

None

pimms.analysis_IO.write_performance(step, eq_string, steps_per_second, overall_moves_per_second, time_elapsed, time_remaining)

Function which writes out the time per step (as taken at some specific step in the simulation) for convenient comparison of simulation efficiency.

Parameters:

• step (int) – The current simulation step number

• eq_string (str) – A string which indicates if the simulation is in equilibrium (E) or production (P) phase.

• steps_per_second (float) – The number of outer master-loop steps per second the simulation is currently running at.

• overall_moves_per_second (float) – The number of individual MC accept/reject moves per second across ALL sub-loops (megamove substeps + TSMMC excursion substeps), i.e. the true MC throughput.

• time_elapsed (str) – The time elapsed in the simulation so far.

• time_remaining (str) – The time remaining in the simulation.

Returns:

No return value, but writes to the file defined in CONFIG.OUTNAME_PERFORMANCE.

Return type:

None

pimms.analysis_IO.write_quench_file(step, temperature, energy)

Append the current step, temperature, and energy to the quench file.

Parameters:

• step (int) – Current simulation step.

• temperature (float) – Current simulation temperature.

• energy (float) – Current total system energy.

Returns:

Nothing is returned; a line is appended to the file defined by CONFIG.QUENCHFILE_NAME.

Return type:

None

pimms.analysis_IO.write_radius_of_gyration(step, RG_list)

Append the instantaneous radius of gyration values for a single step.

Parameters:

• step (int) – Current simulation step.

• RG_list (list of float) – Radius of gyration for each chain (or analysed object) at this step.

Returns:

Nothing is returned; a line is appended to the file defined by CONFIG.OUTNAME_RG.

Return type:

None

pimms.analysis_IO.write_residue_residue_distance(step, R2R_info, all_data)

Append specific residue-residue distance measurements for a single step.

For each monitored residue pair a single line is written containing the step, the two residue indices, and the per-chain distances measured for that pair.

Parameters:

• step (int) – Current simulation step.

• R2R_info (list of sequence of int) – List of residue pairs; each element is an indexable pair where element 0 and element 1 are the two residue indices.

• all_data (list of list of float) – Distance data paired element-wise with R2R_info; each element is a list of distances (e.g. one per chain) for the corresponding pair.

Returns:

Nothing is returned; lines are appended to the file defined by CONFIG.OUTNAME_R2R.

Return type:

None

Raises:

ValueError – If R2R_info and all_data do not have the same length.

pimms.analysis_IO.write_scaling_information(all_nu, all_R0, prefix=False)

Write paired scaling exponent and prefactor values to file.

Writes (overwriting any existing file) a two-column table where each row is a nu (scaling exponent) and R0 (prefactor) pair extracted from the internal-scaling fitting.

Parameters:

• all_nu (list of float) – Scaling exponent values.

• all_R0 (list of float) – Prefactor values, paired element-wise with all_nu.

• prefix (str or bool, optional) – If False (default) the standard output filename is used. If a string is supplied it is prepended to the output basename.

Returns:

Nothing is returned; results are written to the file defined by CONFIG.OUTNAME_SCALING_INFORMATION (optionally prefixed).

Return type:

None

Raises:

ValueError – If all_nu and all_R0 do not have the same length.

class pimms.analysis_structures.DistanceMap(seqlen)

Bases: object

Distance map analysis is an analysis with provides insight into the long-range interaction on a residue-by-residue level - essentially can be considered a contactmap which lacks cutoffs and instead computes the average distance between two residues

get_distance_map()

Return the current system-average distance map.

Returns:

The (seqlen, seqlen) running-mean distance map array (only the upper-right triangle is populated).

Return type:

numpy.ndarray

update_distance_map(dMap)

Fold an instantaneous distance map into the running mean distance map.

Each matrix element is updated as a running average so that the stored matrix always holds the mean over all snapshots seen so far. Only the upper-right triangle is ever populated (the lower triangle stays zero).

Parameters:

dMap (numpy.ndarray) – Square instantaneous distance map with the same shape as the stored matrix.

Return type:

None

Raises:

AnalysisStructureException – If dMap is not a numpy array, or its shape does not match the stored distance map.

write_status(filename='DISTANCE_MAP.dat')

Write the current mean distance map out to file.

Each row is written as comma-separated values, one row per source residue.

Parameters:

filename (str, optional) – Output filename (default 'DISTANCE_MAP.dat'). Overwritten if it already exists.

Return type:

None

class pimms.analysis_structures.InternalScaling(seqlen)

Bases: object

InternalScaling analysis is an analysis with provides insight into the degree of expansion of the chain, but only makes sense in the the context of a full simulation (i.e. an instantaneous Internal Scaling value is not particularly useful).

get_internal_scaling_array()

Return the mean internal scaling values ordered by sequence separation.

Returns:

Mean internal scaling values ordered by increasing gap. The dictionary keys are sorted explicitly because dictionary iteration order is not guaranteed to be numerical.

Return type:

list of float

print_status()

Print the current mean internal scaling profile to stdout.

Return type:

None

update_internal_scaling(IS)

Fold an instantaneous internal scaling profile into the running mean.

Each per-gap value is updated as a running average so that self.internal_scaling always holds the mean over all snapshots seen so far.

Parameters:

IS (dict) – Instantaneous internal scaling values keyed by sequence-separation gap. Must have the same number of entries as the accumulator.

Return type:

None

Raises:

AnalysisStructureException – If IS does not have the same length as the internal accumulator.

write_status(filename='INTSCAL.dat')

Write the current mean internal scaling profile to file.

Parameters:

filename (str, optional) – Output filename (default 'INTSCAL.dat'). Overwritten if it already exists.

Return type:

None

class pimms.analysis_structures.InternalScalingSquared(seqlen)

Bases: object

InternalScalingSquared analysis is an analysis with provides insight into the degree of expansion of the chain, but only makes sense in the the context of a full simulation (i.e. an instantaneous Internal Scaling value is not particularly useful).

fit_scaling_exponent()

Fit the polymer scaling exponent and prefactor from the mean profile.

This method for extracting scaling relationships was developed to avoid the bias introduced by the fact that on a log scale, most inter-residue distances occupy the top-right part of the fitting regime; the idea is to shift to approximately evenly spaced points in log space for the linear fit. The first 15 sequence-separation gaps are always discarded, and at most 40 log-spaced points are used for the fit.

Returns:

(nu, R0) where nu is the fitted scaling exponent and R0 the prefactor. Returns (-1, -1) if the chain is too short (fewer than 25 internal-scaling gaps) to fit meaningfully.

Return type:

tuple of float

get_internal_scaling_array()

Return the mean internal scaling squared values ordered by separation.

Returns:

Mean internal scaling squared values ordered by increasing gap. The dictionary keys are sorted explicitly because dictionary iteration order is not guaranteed to be numerical.

Return type:

list of float

print_status()

Print the current mean internal scaling squared profile to stdout.

Return type:

None

update_internal_scaling(IS)

Fold an instantaneous internal scaling profile into the running mean of the squared distances.

Note that IS is expected to contain the instantaneous internal scaling distances; the value accumulated into the running mean is the square of each distance (IS[i] * IS[i]).

Parameters:

IS (dict) – Instantaneous internal scaling distances keyed by sequence-separation gap. Must have the same number of entries as the accumulator.

Return type:

None

Raises:

AnalysisStructureException – If IS does not have the same length as the internal accumulator.

write_status(filename='INTSCAL.dat')

Write the current mean internal scaling squared profile to file.

Parameters:

filename (str, optional) – Output filename (default 'INTSCAL.dat'). Overwritten if it already exists.

Return type:

None

pimms.lattice_analysis_utils.compute_cluster_gross_properties(cluster_position_list)

Determines the volume of each cluster in a list of clusters. NOTE that positions here MUST have been corrected for PBC effects as the ConvexHull algorithm will calculate ASSUMES a single non-periodic image. This means that when you have clusters that wrap around the PBC the convex hull algorithm is calculating a single instance (i.e. using the boundaries as edges) so take care when extrapolating cluster volume for such system spanning clusters.

Using these positions and the Complex Hull algorithm we compute the volume, area and density of the cluster.

Parameters:

cluster_position_list (list) – list of np.ndarrays where dimensions of the np.ndarray reflect dimensions of the lattice.

Returns:

Returns a list of lists, where each sublist contains three elements that reflect the cluster gross properties. Sublist indices match indices for cluster_position_list indices

[0] - volume [1] - surface area [2] - density

Return type:

list of lists

pimms.lattice_analysis_utils.compute_cluster_radial_density_profile(cluster_position_list, dimensions, minimum_cluster_size_in_beads=None)

Compute the radial density profile of each cluster about its center of mass.

For each cluster the density at “shell k” is the fraction of the lattice sites at Chebyshev (max-norm) distance k from the cluster centre of mass that are occupied by a bead - i.e. (beads at distance k) / (sites in shell k). The profile runs outward from the COM until every bead has been placed in a shell (or the box half-extent is reached), and short profiles are zero-padded to a common length.

This is computed directly by binning each bead’s Chebyshev distance from the COM (an O(num_beads) np.bincount), rather than scanning every site of every concentric shell (which was O(offset_max ** n_dim) and dominated the cost for large clusters). It also fixes an off-by-one in the previous ring-scan, which additionally emitted a spurious shell at offset_max+1 (whose extent spills outside the box); profiles are now capped at offset_max entries as intended.

Parameters:

• cluster_position_list (list of numpy.ndarray) – List of clusters; each entry is an array of lattice positions (each position being a 2- or 3-element coordinate) for that cluster. Positions are expected to be single-image (PBC-corrected) per cluster.

• dimensions (list of int) – Defines the box size in 2 or 3 dimensions (length 2 or 3).

• minimum_cluster_size_in_beads (int or None, optional) – If supplied, clusters with fewer beads than this threshold are skipped (no profile is emitted for them). Default is None (no filtering).

Returns:

One radial density profile per (non-skipped) cluster; each profile is a list of occupied-site fractions as a function of Chebyshev distance from the cluster center of mass, zero-padded to a uniform length.

Return type:

list of list of float

pimms.lattice_analysis_utils.correct_LR_cluster_positions_to_single_image(cluster_position_list, dimensions)

Function which takes a list of cluster positions (i.e. a list of lists, where each sublist is a list of positions associated with the residues in a specific cluster) and for EACH CLUSTER re-configures the cluster position so the cluster is in its own single periodic image

Identical to correct_cluster_positions_to_single_image() but uses a space_threshold of 2, appropriate for long-range (LR) clusters whose members may be separated by more than one lattice site.

Parameters:

• cluster_position_list (list) – A list of lists; each sublist is a list of cluster positions (where each position is itself a 2- or 3-element list).

• dimensions (list) – A list of 2 or 3 elements defining the X/Y or X/Y/Z box dimensions.

Returns:

List of the same length as cluster_position_list where each entry is the cluster’s positions re-expressed in a single (non-periodic) image, using a space_threshold of 2.

Return type:

list

pimms.lattice_analysis_utils.correct_cluster_positions_to_single_image(cluster_position_list, dimensions)

Function which takes a list of cluster positions (i.e. a list of lists, where each sublist is a list of positions associated with the residues in a specific cluster) and for EACH CLUSTER re-configures the cluster position so the cluster is in its own single periodic image.

Parameters:

• cluster_postion_list (list) – A list of lists, each sublist is a list of cluster positions (where, in fact, each position is also itself a list of 2 or 3 positions

• dimensions (list) – A list of 2- or 3- elements that defines the X/Y or X/Y/Z positions

Returns:

List of the same length as cluster_position_list where each entry is the cluster’s positions re-expressed in a single (non-periodic) image, as returned by cluster_utils.convert_positions_to_single_image_snakesearch with a space_threshold of 1.

Return type:

list

pimms.lattice_analysis_utils.extract_cluster_polymeric_properties(cluster_position_list, dimensions=False)

Function which takes a list of cluster positions (i.e. a list of lists, where each sublist is a list of positions associated with the residues in a specific cluster) and returns a list of lists of the same length where each sublist in the return list contains the polymeric properties of the actual cluster.

Parameters:

• cluster_position_list (list of list of positions) – List where each sublist is a list of positions; each sublist is its own cluster. NOTE that each cluster should exist within its own single-image convention, so that for each cluster the properties can be computed naively over those positions without any further PBC handling.

• dimensions (list of int or bool, optional) – Defines the dimensions of the lattice the positions sit on. If PBC correction has already been performed this can be left as False (the default), in which case a per-cluster bounding box is computed dynamically (+10 beyond the largest value in each dimension).

Returns:

List of the same length as cluster_position_list where each entry is the [rg, asph] polymeric properties of the corresponding cluster (empty list if no clusters are supplied).

Return type:

list of list of float

pimms.lattice_analysis_utils.extract_positions_from_clusters(cluster_list, chainDict)

Function which takes a list of clusters (i.e. a list of lists, where each sublist is a list of chainIDs in a specific cluster) and returns a list of lists of the same length where each sublist in the return list contains the positon of all residues in the cluster

Parameters:

• cluster_list (list of list of int) – List of clusters, where each sublist is a list of chainIDs in that cluster.

• chainDict (dict) – Dictionary mapping each chainID to its chain object (each chain object must expose get_ordered_positions()).

Returns:

List of the same length as cluster_list where each sublist contains the ordered positions of all residues belonging to that cluster.

Return type:

list of list

pimms.lattice_analysis_utils.get_LR_cluster_distribution(latticeObject)

Returns a list of lists, where each sublist contains the chainIDs associated with a cluster. Cluster sublists are ordered from largest cluster to smallest. LR clusters are defined as clusters were interactions are through short-range OR long-range interactions

Parameters:

latticeObject (Lattice) – The lattice object. Its grid (the lattice grid) and chains (mapping of chainIDs to chain objects) attributes are used, along with long-range interaction information, to build the long-range clusters.

Returns:

List of clusters, where each sublist contains the chainIDs of the chains in that long-range connected component. Clusters are ordered from largest to smallest.

Return type:

list of list of int

pimms.lattice_analysis_utils.get_cluster_distribution(lattice_grid, chainDict)

Returns a list of lists, where each sublist contains the chainIDs associated with a cluster. Cluster sublists are ordered from largest cluster to smallest.

This is a computationally expensive algorithm that probably could be ported into Cython at some point…

Parameters:

• lattice_grid (np.array (2D or 3D)) – Standard lattice grid

• chainDict (dict) – Standard dicionary mapping chainIDs to chain objects.

Returns:

List of clusters, where each sublist contains the chainIDs of the chains in that connected component. Clusters are ordered from largest to smallest.

Return type:

list of list of int

pimms.lattice_analysis_utils.get_eigenvalues_of_the_T_matrix(positions, dimensions, pbc_correction=True)

Compute the eigenvalues and eigenvectors of the gyration (T) tensor.

Builds the gyration tensor from the supplied positions relative to their (optionally PBC-corrected) center of mass, then diagonalizes it. The eigenvalues are the principal components used downstream to compute the radius of gyration and asphericity.

Parameters:

• positions (list of positions) – A list of positions, where each position is a 2-length or 3-length list specifying X/Y/[Z] coordinate positions on the lattice.

• dimensions (list of int) – Defines the box size in 2 or 3 dimensions (length 2 or 3).

• pbc_correction (bool, optional) – If True (default), each position is PBC-corrected relative to the center of mass before contributing to the gyration tensor.

Returns:

(EIG, norm) where EIG is the array of eigenvalues of the gyration tensor and norm is the matrix of corresponding eigenvectors (as returned by numpy.linalg.eig).

Return type:

tuple

pimms.lattice_analysis_utils.get_inter_position_distance(P1, P2, dimensions, pbc_correction=True)

Returns the distance between two positions on the lattice (in real space) accounting for periodic boundary conditions.

Routine optimized for a single distance (i.e. doesn’t perform any of the setup/teardown used for vectorized implementations which are important when a set of positions are being compared)

Parameters:

• P1 ([list of ints]) – A position list (e.g. a list of integers specifying the X/Y or X/Y/Z coordinates of a position)

• P2 ([list of ints]) – A position list (e.g. a list of integers specifying the X/Y or X/Y/Z coordinates of a position)

• dimensions ([list of ints, 2 or 3 in length]) – Defines the box size in 2 or 3 dimensions

• pbc_correction (bool) – Flag which if set to true means a PBC correction is applied. Default is True.

Returns:

The (optionally PBC-corrected) Euclidean distance between P1 and P2 in real space.

Return type:

float

pimms.lattice_analysis_utils.get_inter_position_distances(P1s, P2s, dimensions, pbc_correction=True)

Returns the list of distances between lists of two positions on the lattice (in real space) accounting for periodic boundary conditions.

Optimized for multiple values - vectorizes calculations.

Parameters:

• P1s (list of positions) – A list of positions, where each position is a 2-length or 3-length list specifying X/Y/[Z] coordinate positions on the lattice.

• P2s (list of positions) – A list of positions, where each position is a 2-length or 3-length list specifying X/Y/[Z] coordinate positions on the lattice.

• dimensions (list of int) – Defines the box size in 2 or 3 dimensions (length 2 or 3).

• pbc_correction (bool, optional) – If True (default), the minimum-image PBC correction is applied to each per-dimension separation before computing distances.

Returns:

1D array of (optionally PBC-corrected) Euclidean distances, one per position pair.

Return type:

numpy.ndarray

Raises:

AnalysisRoutineException – If P1s and P2s do not have the same length.

pimms.lattice_analysis_utils.get_polymeric_properties(positions, dimensions, pbc_correction=True)

Returns a list of polymeric properties calculated over the set of positions

[0] - radius of gyration [1] - asphericity

Rg is defined as

sqrt(dfrac{1}{N}sum_{k=1}^N(r_k-r_{mean})^2)

Where N = number of residues r_{mean} = mean residue position (Center of Mass)

Parameters:

• positions (list of positions) – A list of positions, where each position is a 2-length or 3-length list specifying X/Y/[Z] coordinate positions on the lattice.

• dimensions (list of int) – Defines the box size in 2 or 3 dimensions (length 2 or 3).

• pbc_correction (bool, optional) – Defines whether to perform PBC correction here (default True). For certain types of analysis (notably cluster analysis) the PBC correction is dealt with by the algorithms that construct the cluster, such that performing it again here is redundant (and generally not possible, as the snakesearch algorithm re-positions the cluster in terms of non-periodic space).

Returns:

A two-element list [rg, asph] where rg is the radius of gyration and asph is the asphericity (acylindricity in 2D), both derived from the gyration-tensor eigenvalues. Degenerate cases where rg ~ 0 return an asphericity of 0.0.

Return type:

list of float