Forces

Overview

Bases: object

Base class from which all other forces inherit. Don’t use this class directly, instead use Bond, Angle, Pair, and Dihedral.

Warning

This class should not be instantiated directly by users. It can be used for isinstance or issubclass.

Note

Forces in MSIBI can either be held constant (i.e., fixed) or optimized (i.e., mutable).

Only one type of force can be optimized at a time. For example, you can optimize multiple angle potentials during one optimization run, but you cannot optimize a Pair and an Angle potential in the same optimization run.

Several of the methods in this class are only applicable in the case a force is being optimized.

Parameters:

name (str) – The name of the type in the Force. Must match the names found in the State’s trajectory file.
optimize (bool) – Set to True if this force is to be mutable and optimized. Set to False if this force is to be held constant while other forces are optimized.
nbins (int, optional) – This must be a positive integer if this force is being optimized. nbins is used for setting the potenials independent varible (x) range and step size (dx). It is also used in determining the bin size of the target and query distributions. If this force is not being optimied, leave this as None.
smoothing_window (int, optional) – The window size used in SciPy’s savgol_filter method.
smoothing_order (int, optional) – The smoothing order used in SciPy’s savgol_filter method.
correction_fit_window (int, optional) – The window size (number of data points) to use when fitting the iterative potential to head and tail correction forms. This is only used when the Force is set to be optimized.
maxfev (int, optional) – Sets the maximum number of attemps when using scipy.curve_fit to apply head and tail corrections. Use larger values to improve likelihood of successful corrections, but this may slow performance.
correction_form (Callable, optional) – The type of correciton form to apply to the potential. This is only used when the Force is set to be optimized. Bonded forces (Bond, Angle, Dihedral) apply this correction to both the head and tail of the potential. Non-bonded forces (Pair) only apply this correction to the head of the potential.

property potential: ndarray: The potential energy values V(x).

property force: ndarray: The force values F(x).

property smoothing_window: int: Window size used in smoothing the distributions and potentials.

property smoothing_order: int: The order used in Savitzky Golay filter.

property nbins: int: The number of bins used in distributions and x-range of the potential.

smooth_potential() → None: Smooth and overwrite the current potential.

Note

This uses a Savitzky Golay smoothing algorithm where the window size and order parameters are set by msibi.forces.Force.smoothing_window() and msibi.forces.Force.smoothing_order().

Both of these can be changed using their respective setters.

save_potential(file_path: str) → None

Save the x-range, potential and force to a .csv file.

Note

This method uses pandas.DataFrame.to_csv and saves the data in with column labels of “x”, “potential”, and “force”.

If you want to smooth the final potential, use msibi.forces.Force.smooth_potential() before calling this method.

Parameters:: file_path (str) – File path and name to save table potential to.

save_potential_history(file_path: str) → None

Save the potential history of the force to a .npy file.

Parameters:: file_path (str) – File path and name to save table potential history to.

save_state_data(state: State, file_path: str) → None

Save the distribution data of a state as a a dictionary to a .npz file.

Note

This saves the state points target distribution, current distribution distribution history and f-fit scores for the corresponding potential to the .npz file.

Parameters:

state (msibi.state.State, required) – The state to use in finding the target distribution.
file_path (str, required) – File path and name to save the .npz file.

target_distribution(state: State) → ndarray

The target structural distribution corresponding to this force.

Parameters:: state (msibi.state.State) – The state point to use in finding the target distribution.

plot_target_distribution(state: State, file_path: str = None) → None

Plot the target distribution corresponding to this force and state point.

Note

Use this to see how the shape of the target distribution is affected by your choices for nbins, smoothing window, and smoothing order.

Parameters:

state (msibi.state.State) – The state to use in finding the target distribution.
file_path (str, optional) – If given, the plot will be saved to this location.

plot_fit_scores(state: State, file_path: str = None) → None

Plot the evolution of the distribution matching fit scores.

Parameters:

state (msibi.state.State) – The state to use in finding the target distribution.
file_path (str, optional) – If given, the plot will be saved to this location.

plot_potential(file_path: str | None = None, xlim: tuple | None = None, ylim: tuple | None = None) → None

Plot the currently optimized potential energy.

Parameters:

file_path (str, optional) – If given, the plot will be saved to this location.
xlim (tuple, optional) – If given, sets the limits for the x-range of the plot. If not given, uses the entire x-range used in optimization.
ylim (tuple, optional) – If given, sets the limits for the y-range of the plot. If not given, uses the entire y-range of the learned potential.

plot_potential_history(file_path: str | None = None, xlim: tuple | None = None, ylim: tuple | None = None) → None

Plot the history of the optimized potential energy.

Parameters:

file_path (str, optional) – If given, the plot will be saved to this location.
xlim (tuple, optional) – If given, sets the limits for the x-range of the plot. If not given, uses the entire x-range used in optimization.
ylim (tuple, optional) – If given, sets the limits for the y-range of the plot. If not given, uses the entire y-range of the learned potential.

plot_distribution_comparison(state: State, file_path: str | None = None) → None

Plot the target distribution and most recent query distribution.

Parameters:

state (msibi.state.State) – The state to use in finding the target distribution.
file_path (str, optional) – If given, the plot will be saved to this location.

distribution_history(state: State) → ndarray

Get the complete query distribution history for a given state.

Parameters:: state (msibi.state.State) – The state point to use for calculating the distribution.

set_target_distribution(state: State, array: ndarray) → None

Manually set the target distribution for a given state.

Parameters:

state (msibi.state.State) – The state point used in finding the distribution.
array (np.ndarray) – The 2D array representing the target distribution for this state point.

current_distribution(state: State) → ndarray

Distrubution of the most recent iteration for a given state point.

Parameters:: state (msibi.state.State) – The state point used for calculating the distribution.

distribution_fit(state: State) → float

Get the fit score from the most recent query simulation.

Parameters:: state (msibi.state.State) – The state point used for calculating the distribution.

Set a potential based on the following function:

\(V(x) = k4(x-x_{0})^{4} + k3(x-x_{0})^{3} + k2(x-x_{0})^{2}\)

Note

Using this method will create a table potential V(x) over the range x_min - x_max.

This is useful for easily setting initial guess potentials for bond, angle and dihedral forces to be optimized. Also see msibi.forces.Force.set_from_file().

See msibi.forces.Pair.set_lj() for setting an initial guess potential for a non-bonded pair force to be optimized.

Parameters:

x0 (float) – The paraters used in the V(x) function described above.
k2 (float) – The paraters used in the V(x) function described above.
k3 (float) – The paraters used in the V(x) function described above.
k4 (float) – The paraters used in the V(x) function described above.
x_min (float) – The lower bound of the potential range.
x_max (float) – The upper bound of the potential range.

set_from_file(file_path: str) → None

Set a potential from a .csv file.

Warning

This uses pandas.DataFrame.read_csv and expects column names of “x”, “potential”, and “force”.

Tip

Use this potential setter to set a potential from a previous MSIBI run. For example, use the final potential files from a bond-optimization IBI run to set a static coarse-grained bond potential while you perform IBI runs on angle and/or pair potentials.

Also see: msibi.forces.Force.save_potential()

Parameters:: file_path (str) – The full path to the table potential csv file.

class msibi.forces.Bond(type1: str, type2: str, optimize: bool, nbins: int | None = None, smoothing_window: int | None = 15, smoothing_order: int | None = 2, correction_fit_window: int | None = 10, maxfev: int | None = 1000, correction_form: ~typing.Callable = <function harmonic>)

Bases: Force

Creates a bond stretching Force used in query simulations.

Note

The bond type is sorted so that type1 and type2 are listed in alphabetical order, and must match the bond types found in the state’s target GSD file bond types.

For example: Bond(type1="B", type2="A") will have Bond.name = "A-B"

Parameters:

type1 (str) – Name of the first particle type in the bond. This must match the types found in the state’s target GSD file.
type2 (str) – Name of the second particle type in the bond. This must match the types found in the state’s target GSD file.
optimize (bool) – Set to True if this force is to be mutable and optimized. Set to False if this force is to be held constant while other forces are optimized.
nbins (int, optional) – This must be a positive integer if this force is being optimized. nbins is used for setting the potenials independent varible (x) range and step size (dx). It is also used in determining the bin size of the target and query distributions.
smoothing_window (int, optional default 15) – The window size used in SciPy’s savgol_filter method.
smoothing_order (int, optional default 2) – The smoothing order used in SciPy’s savgol_filter method.
correction_fit_window (int, optional default 10) – The window size (number of data points) to use when fitting the iterative potential to head and tail correction forms. This is only used when the Force is set to be optimized.
maxfev (int, optional default 3000) – Sets the maximum number of attemps when using scipy.curve_fit to apply head and tail corrections. Use larger values to improve likelihood of successful corrections, but this may slow performance.
correction_form (Callable, default=msibi.utils.corrections.harmonic) – The type of correciton form to apply to the potential. This is only used when the Force is set to be optimized. This correction is applied to both the head and tail of the potential.

set_state_params(state: State, optimize_against: bool) → None

Set Bond-State specific parameters.

Note

Use this method to override default values for this force, which are used for all states included the MSIBI optimization.

One use case is ignore a certain State for a Force that is being optimized.

Parameters:

state (msibi.state.State) – The state point being used to define Force-state specific parameters.
optimize_against (bool) – Turns this state off (False) or on (True), for this Force. This overrides Force.optimize only for this state point.

set_harmonic(r0: float | int, k: float | int) → None

Set a fixed harmonic bond potential.

Warning

Using this method is not compatible with Force objects that are set to be optimized during MSIBI.

Note

For more information on harmonic bond potentials, refer to the HOOMD-blue harmonic bond documentation.

Parameters:

r0 ((Union[float, int])) – Equilibrium bond length [length]
k ((Union[float, int])) – Spring constant [energy]

class msibi.forces.Angle(type1: str, type2: str, type3: str, optimize: bool, nbins: int | None = None, smoothing_window: int | None = 15, smoothing_order: int | None = 2, correction_fit_window: int | None = 10, maxfev: int | None = 1000, correction_form: ~typing.Callable = <function harmonic>)

Bases: Force

Creates a bond angle Force used in query simulations.

Note

The angle type is formed in the order of type1-type2-type3 and must match the same order in the target GSD file angle types.

Parameters:

type1 (str) – Name of the first particle type in the angle. This must match the types found in the state’s target GSD file.
type2 (str) – Name of the second particle type in the angle. This must match the types found in the state’s target GSD file.
type3 (str) – Name of the third particle type in the angle. This must match the types found in the state’s target GSD file.
optimize (bool) – Set to True if this force is to be mutable and optimized. Set to False if this force is to be held constant while other forces are optimized.
nbins (int, otional) – This must be a positive integer if this force is being optimized. nbins is used for setting the potenials independent varible (x) range and step size (dx). It is also used in determining the bin size of the target and query distributions.
smoothing_window (int, optional default 15) – The window size used in SciPy’s savgol_filter method.
smoothing_order (int, optional default 2) – The smoothing order used in SciPy’s savgol_filter method.
correction_fit_window (int, optional default 10) – The window size (number of data points) to use when fitting the iterative potential to head and tail correction forms. This is only used when the Force is set to be optimized.
maxfev (int, optional) – Sets the maximum number of attemps when using scipy.curve_fit to apply head and tail corrections. Use larger values to improve likelihood of successful corrections, but this may slow performance.
correction_form (Callable, default=msibi.utils.corrections.harmonic) – The type of correciton form to apply to the potential. This is only used when the Force is set to be optimized. This correction is applied to both the head and tail of the potential.

set_state_params(state: State, optimize_against: bool) → None

Set Angle-State specific parameters.

Note

Use this method to override default values for this force, which are used for all states included the MSIBI optimization.

One use case is ignore a certain State for a Force that is being optimized.

Parameters:

state (msibi.state.State) – The state point being used to define Force-state specific parameters.
optimize_against (bool) – Turns this state off (False) or on (True), for this Force. This overrides Force.optimize only for this state point.

set_harmonic(t0: float | int, k: float | int) → None

Set a fixed harmonic angle potential.

Warning

Using this method is not compatible with Force objects that are set to be optimized during MSIBI.

Note

For more information on harmonic angle potentials, refer to the HOOMD-blue harmonic angle documentation.

Parameters:

t0 ((Union[float, int])) – Equilibrium bond angle [radians]
k ((Union[float, int])) – Spring constant [energy]

class msibi.forces.Pair(type1: str, type2: str, optimize: bool, nbins: int | None = None, r_cut: float | int | None = None, r_switch: float | int | None = None, exclude_bond_depth: int = 0, exclude_all_bonded: bool = False, smoothing_window: int | None = 11, smoothing_order: int | None = 2, correction_fit_window: int | None = 8, maxfev: int | None = 1000, head_correction_form: ~typing.Callable = <function exponential>)

Bases: Force

Creates a non-bonded pair Force used in query simulations.

Note

The pair type is sorted so that type1 and type2 are listed in alphabetical order, and must match the pair types found in the state’s target GSD file bond types.

For example: Pair(type1="B", type2="A") will have Pair.name = "A-B"

Parameters:

type1 (str) – Name of the first particle type in the pair. This must match the types found in the state’s target GSD file.
type2 (str) – Name of the second particle type in the pair. This must match the types found in the state’s target GSD file.
optimize (bool) – Set to True if this force is to be mutable and optimized. Set to False if this force is to be held constant while other forces are optimized.
r_cut ((Union[float, int])) – Sets the cutoff distance used in Hoomd’s neighborlist. This is also the maximum radius value used in the pair table potential.
r_switch ((Union[float, int])) – Sets the radius value where the potential is corrected to begin approach zero smoothly.
nbins (int, optional) – This must be a positive integer if this force is being optimized. nbins is used for setting the potenials independent varible (x) range and step size (dx). It is also used in determining the bin size of the target and query distributions.
exclude_bond_depth (int, optional (default 0)) – Excludes all pairs within a depth (distance on a bond graph) from the RDF calculation. This must be a positive integer.
exclude_all_bonded (bool, optional (default False)) – Excludes all pairs belonging to the same molecule from the RDF calculation. If this is used, then exclude_bond_depth should not be used.
smoothing_window (int, optional default 11) – The window size used in SciPy’s savgol_filter method.
smoothing_order (int, optional default 2) – The smoothing order used in SciPy’s savgol_filter method.
correction_fit_window (int, optional default 8) – The window size (number of data points) to use when fitting the iterative potential to head and tail correction forms. This is only used when the Force is set to be optimized.
maxfev (int, optional) – Sets the maximum number of attemps when using scipy.curve_fit to apply head and tail corrections. Use larger values to improve likelihood of successful corrections, but this may slow performance.
correction_form (Callable, default=msibi.utils.corrections.harmonic) – The type of correciton form to apply to the potential. This is only used when the Force is set to be optimized. This correction is only applied to both the head of the potential.

set_state_params(state: State, optimize_against: bool, exclude_bond_depth: int, exclude_all_bonded: bool) → None

Set Pair-State specific parameters.

Note

Use this method to override default values for this force, which are used for all states included the MSIBI optimization.

One use case is ignore a certain State for a Force that is being optimized.

Parameters:

state (msibi.state.State) – The state point being used to define Force-state specific parameters.
optimize_against (bool) – Turns this state off (False) or on (True), for this Force. This overrides Force.optimize only for this state point.
exclude_bond_depth (int) – Sets a bond-depth pair exclusion rule that applies only to calculating RDFs for this state.
exclude_all_bonded (bool) – Excludes all intra-molecular pairs in the RDF calculation

set_lj(r_min: float | int, r_cut: float | int, epsilon: float | int, sigma: float | int) → None

Set a 12-6 Lennard Jones table pair potential used in query simulations.

Note

This creates a table potential from the LJ 12-6 function with the given parameters. Use this to create an initial guess when optimizing a Pair force. It can still be used to set a static potential.

Parameters:

epsilon ((Union[float, int])) – Sets the dept hof the potential energy well.
sigma ((Union[float, int])) – Sets the particle size.
r_cut ((Union[float, int])) – Maximum distance used to calculate neighbors.

class msibi.forces.Dihedral(type1: str, type2: str, type3: str, type4: str, optimize: bool, nbins: int | None = None, smoothing_window: int | None = 11, smoothing_order: int | None = 2, correction_fit_window: int | None = 10, maxfev: int | None = 1000, correction_form: ~typing.Callable = <function harmonic>)

Bases: Force

Creates a dihedral Force used in query simulations.

Note

The dihedral type is formed in the order of type1-type2-type3-type4 and must match the same order in the target GSD file dihedral types.

Parameters:

type1 (str) – Name of the first particle type in the dihedral. This must match the types found in the state’s target GSD file.
type2 (str) – Name of the second particle type in the dihedral. This must match the types found in the state’s target GSD file.
type3 (str) – Name of the third particle type in the dihedral. This must match the types found in the state’s target GSD file.
type4 (str) – Name of the fourth particle type in the dihedral. This must match the types found in the state’s target GSD file.
optimize (bool) – Set to True if this force is to be mutable and optimized. Set to False if this force is to be held constant while other forces are optimized.
nbins (int, optional) – This must be a positive integer if this force is being optimized. nbins is used for setting the potenials independent varible (x) range and step size (dx). It is also used in determining the bin size of the target and query distributions.
smoothing_window (int, optional default 11) – The window size used in SciPy’s savgol_filter method.
smoothing_order (int, optional default 2) – The smoothing order used in SciPy’s savgol_filter method.
correction_fit_window (int, optional default 10) – The window size (number of data points) to use when fitting the iterative potential to head and tail correction forms. This is only used when the Force is set to be optimized.
maxfev (int, optional) – Sets the maximum number of attemps when using scipy.curve_fit to apply head and tail corrections. Use larger values to improve likelihood of successful corrections, but this may slow performance.
correction_form (Callable, default=msibi.utils.corrections.harmonic) – The type of correciton form to apply to the potential. This is only used when the Force is set to be optimized. This correction is applied to both the head and tail of the potential.

set_state_params(state: State, optimize_against: bool) → None

Set Dihedral-State specific parameters.

Note

Use this method to override default values for this force, which are used for all states included the MSIBI optimization.

One use case is ignore a certain State for a Force that is being optimized.

Parameters:

state (msibi.state.State) – The state point being used to define Force-state specific parameters.
optimize_against (bool) – Turns this state off (False) or on (True), for this Force. This overrides Force.optimize only for this state point.

set_periodic(phi0: float | int, k: float | int, d: int, n: int) → None

Set a fixed periodic dihedral potential.

Warning

Using this method is not compatible with Force objects that are set to be optimized during MSIBI.

Note

For more information on periodic dihedral potentials, refer to the HOOMD-blue dihedral periodic documentation.

Parameters:

phi0 (float) – Phase shift [radians]
k (float) – Spring constant [energy]
d (int) – Sign factor
n (int) – Angle multiplicity