Options

There are some options to control what exactly is going to be computed to obtain the MDDF. These options can be defined by the user and passed to the mddf function, using, for example:

options = Options(lastframe=1000, cutoff=12.0)
results = mddf(trajectory_file, solute, solvent, options)

Frame ranges and histogram properties

These are common options that the regular user might want to set in their calculation.

firstframe: Integer, first frame of the trajectory to be considered.

lastframe: Integer, last frame of the trajectory to be considered.

stride: Integer, consider every stride frames, that is, if stride=5 only one in five frames will be considered.

binstep: Real, length of the bin step of the histograms, default = 0.02 Angstroms.

cutoff: Real, the maximum distance to the solute considered in the construction of the histograms. The bulk density of the solvent is estimated from the molecules beyond the cutoff, by exclusion: from the number of solvent molecules beyond the cutoff and from the volume of this region, computed by numerical integration. Set this distance to a value beyond which the distribution function is converged to 1 for practical purposes. Default: 10 Angstroms, but setting it according to the system is highly encouraged. The dependence of the density of the solvent on the distance to the solute, beyond the correlation length, can be inspected with the reference_density function.

Lower level options

These will probably never be set by the user, unless if dealing with some special system (large very large, or very low density system).

irefatom: Integer, index of the reference atom in the solvent molecule used to compute the shell volumes and domain volumes in the Monte-Carlo volume estimates. The final rdf data is reported for this atom as well. By default, we choose the atom which is closer to the center of coordinates of the molecule, but any choice should be fine.

n_random_samples: Integer, the number of random (ideal-gas) samples computed for each frame, used to compute the shell volumes and the random MDDF counts that normalize the distributions. The random samples are generated as follows:

• For cross-correlations (the solute and solvent are different molecules) with more than one solute molecule, all solvent molecules are placed with random positions and orientations (keeping their internal coordinates), and the minimum distances of the random solvent molecules to every solute molecule are computed, as for the actual frame. Thus, each sample provides one reference for each solute molecule. The default is 1 sample per frame.

• For cross-correlations with a single solute molecule (for example, a protein), all solvent molecules are placed with random positions and orientations (with internal coordinates copied from solvent molecules in the bulk), and the minimum distances of the random solvent molecules to the solute molecule are computed. The default is 10 samples per frame.

• For self-correlations (the solute and solvent are the same molecules), all molecules are placed with random positions and orientations (keeping their internal coordinates), and the minimum distances of each molecule to all other molecules are computed, exactly as for the actual frame. Thus, each sample is a random configuration of the system, which provides one reference molecule for each molecule of the system. The default is 1 sample per frame.

The number of random samples per frame should be large enough such that the random counts are sampled at least as well as the actual minimum-distance counts. The defaults satisfy this condition. Increase n_random_samples if you have short trajectories and want to obtain reproducible results for them.

seed: Seed for random number generator. If -1, the seed will be generated from the entropy of the system. If your results are dependent on the seed, it is probable that you do not have enough sampling. Mostly used for testing purposes. Two runs are only identical if ran with the same seed and in serial mode.

StableRNG (::Bool), defaults to false. Use a stable random number generator from the StableRNGs package, to produce identical runs on different architectures and Julia versions. Only used for testing.

nthreads: How many threads to use. By default, it will be the number of physical cores of the computer.

lcell: Integer, the cell length of the linked-cell method (actually the cell length is cutoff/lcell). Default: 1.

GC: Bool, force garbage collection, to avoid memory overflow. Default: true. That this might be required is probably a result of something that can vastly improved in memory management. This may slow down parallel runs significantly if the GC runs too often.

GC_threshold: Float64, minimum fraction of the total memory of the system required to force a GC run. That is, if GC_threshold=0.1, which is the default, every time the free memory becomes less or equal to 10% of the total memory available, a GC run occurs.

Frame statistical reweighing

Most times the weights of each frame of the trajectory are the same, resulting from some standard MD simulation. If, for some reason, the frames have different statistical weights, the weights can be passed to the as an optional parameter frame_weights.

For example:

julia> results = mddf(trajectory_file, solute, solvent, options; frame_weights=[0.0, 1.0, 2.0])

The code above will assign a larger weight to the third frame of the trajectory. These weights are relative (meaning that [0.0, 1.0, 2.0] would produce the same result). What will happen under the hood is that the distance counts of the frames will be multiplied by each frame weight, and normalized for the sum of the weights.

Important: The length of the frame_weights vector must be at least equal to the number of the last frame read from the trajectory. That is, if lastframe is not set, and all the frames will be read, the length of frame_weights must be equal to the length of the trajectory (the stride parameter will skip the information both of the frames and its weights). If lastframe is set, then the length of frame_weights must be at least lastframe (it can be greater, and further values will be ignored). Importantly, the indices of the elements in frame_weights are assumed to correspond to the indices of the frames in the original trajectory file.

Compute coordination number only

For some systems, it may be impossible, or to expensive, to compute the normalization of the minimum-distance distribution function. Nevertheless, the coordination number may still be an interesting information to be retrieved from the simulations. To run the computation to compute coordination numbers only, do:

julia> results = mddf(trajectory_file, solute, solvent, options; coordination_number_only = true)

Options(;
    firstframe::Int = 1,
    lastframe::Int = -1,
    stride::Int = 1,
    irefatom::Int = -1,
    n_random_samples::Union{Nothing,Int} = nothing,
    binstep::Float64 = 0.02,
    cutoff::Union{Nothing,Real} = nothing,
    lcell::Int = 1,
    GC::Bool = true,
    GC_threshold::Float64 = 0.3,
    seed::Int = 321,
    StableRNG::Bool = false,
    nthreads::Int = 0,
    silent::Bool = false
)

Loading trajectories

To initialize a trajectory file for computation before calling the mddf or coordination_number functions, use the command

trajectory = Trajectory("trajectory.xtc",solute,solvent)

where solute and solvent are defined with the AtomSelection function described before. This function opens the stream for reading frames, which are read once a time when the coordinates are required for computing the MDDF.

This object can be used to feed mddf and coordination_number with:

mddf(trajectory, options)
coordination_number(trajectory, options)

The Trajectory function uses Chemfiles in background, and thus the most common trajectory formats are supported, as the ones produced with NAMD, Gromacs, LAMMPS, Amber, etc.

trajectory = Trajectory("trajectory.xtc",solute,solvent,format="xtc")

Trajectory(filename::AbstractString, solute::AtomSelection, solvent::AtomSelection; format::AbstractString = "", chemfiles = false)

Stop a computation

To stop a running calculation without losing the Julia REPL state, add a file named

stop_complexmixtures

in the same directory from where the computation was launched. This will interrupt the computations in each parallel task as soon as the file is detected. The partial results will be returned and are, thus, invalid. Remove the file to restart a calculation. The content of the file is not important, thus it can be created by the touch stop_complexmixtures linux command.