pyMFI package

Submodules

pyMFI.MFI module

pyMFI.MFI.FFT_intg_2D(FX, FY, min_grid=array([- 3.14159265, - 3.14159265]), max_grid=array([3.14159265, 3.14159265]), periodic=array([0, 0]))

2D integration of force gradient (FX, FY) to find FES using Fast Fourier Transform.

Parameters

• FX (array of size (nbins[1], nbins[0])) – CV1 component of the Mean Force.

• FY (array of size (nbins[1], nbins[0])) – CV1 component of the Mean Force.

• min_grid (array, optional) – Lower bound of the simulation domain. Defaults to np.array((-np.pi, -np.pi)).

• min_grid – Upper bound of the simulation domain. Defaults to np.array((np.pi, np.pi)).

• periodic (list or array of shape (2,))) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

Returns

[X, Y, fes]

X (array of size (nbins[1], nbins[0])): CV1 grid positions

Y (array of size (nbins[1], nbins[0])): CV2 grid positions

fes (array of size (nbins[1], nbins[0])): Free Energy Surface

Return type

list

pyMFI.MFI.MFI_2D(HILLS='HILLS', position_x='position_x', position_y='position_y', bw=array([0.1, 0.1]), kT=1, min_grid=array([- 3.14159265, - 3.14159265]), max_grid=array([3.14159265, 3.14159265]), nbins=array([200, 200]), error_pace=- 1, base_terms=0, window_corners=[], WellTempered=1, nhills=- 1, periodic=array([0, 0]), Ftot_den_limit=1e-10, FES_cutoff=- 1, Ftot_den_cutoff=0.1, non_exploration_penalty=0, use_weighted_st_dev=True, hp_centre_x=0.0, hp_centre_y=0.0, hp_kappa_x=0, hp_kappa_y=0, lw_centre_x=0.0, lw_centre_y=0.0, lw_kappa_x=0, lw_kappa_y=0, uw_centre_x=0.0, uw_centre_y=0.0, uw_kappa_x=0, uw_kappa_y=0, F_static_x=array([[0.]]), F_static_y=array([[0.]]), ref_fes=array([[0.]]))

Compute a time-independent estimate of the Mean Thermodynamic Force, i.e. the free energy gradient in 2D CV spaces.

Parameters

• HILLS (array) – HILLS array. The colums have format: [time [ps], position_x [nm], position_y [nm], MetaD_sigma_x [nm], MetaD_sigma_y [nm], MetaD_height [nm], MetaD_biasfactor]

• position_x (array) – CV1 array. Defaults to “position_x”.

• position_y (array) – CV2 array. Defaults to “position_y”.

• bw (list or array of shape (2,), optional) – Scalar, bandwidth for the construction of the KDE estimate of the biased probability density. First entry is the bandwidth for CV1 and second entry is the bandwidth for CV2. Defaults to np.array((0.1,0.1)).

• kT (float, optional) – Boltzmann constant multiplied with temperature (reduced format, 120K -> 1).

• min_grid (array, optional) – Lower bound of the force domain. Defaults to np.array((-np.pi, -np.pi)).

• max_grid (array, optional) – Upper bound of the force domain. Defaults to np.array((np.pi, np.pi)).

• nbins (array, optional) – number of bins in CV1,CV2. First enrty is the number of bins in CV1 and the second entry is the number of bins in CV2! Defaults to np.array((200,200)).

• error_pace (int, optional) – Pace for the calculation of the on-the-fly measure of global convergence. Defaults to 1, change it to a higher value if FES_cutoff>0 is used.

• base_terms (int or list, optional) – When set to 0, inactive. When activated, “on the fly” variance is calculated as a patch to base (previous) simulation. To activate, put force terms of base simulation ([Ftot_den, Ftot_den2, Ftot_x, Ftot_y, ofv_x, ofv_y]). Defaults to 0.

• window_corners (list, optional) – When set to [], inactive. When activated, error is ALSO calculated for mean force in the window. To activate, put the min and max values of the window ([min_x, max_x, min_y, max_y]). Defaults to [].

• WellTempered (binary, optional) – Is the simulation well tempered? 1 for yes and 0 for no. Defaults to 1.

• nhills (int, optional) – Number of HILLS to analyse, -1 for the entire HILLS array. Defaults to -1, i.e. the entire dataset.

• periodic (list or array of shape (2,), optional) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1. Defaults to np.array((0,0)).

• Ftot_den_limit (float, optional) – Truncation of the probability density for numerical reasons, to avaiod devisions by zero (or suare root of negative numbers). If the probability density (Ftot_den) of some CV region is lover than the Ftot_den_limit, it will be set to zero. Default is set to 1E-10.

• FES_cutoff (float, optional) – Cutoff applied to error calculation for FES values over the FES_cutoff. If the cutoff applies, the error will be set to zero, otherwise the error will stay the same. Defaults to 0. When FES_cutoff <= 0, no cufoff is applied. Use with care, computing the fes in the loop renders the calculation slow.

• Ftot_den_cutoff (float, optional) – Cutoff applied to error calculation for probability density (Ftot_den) values below the Ftot_den_cutoff. If the cutoff applies, the error will be set to zero, otherwise the error will stay the same. Defaults to 0.1. When Ftot_den_cutoff <= 0, no cufoff is applied.

• non_exploration_penalty (float, optional) – Turns zero-value error to the non_exploration_penalty value. This should be used in combination with the cutoff. If some part of CV space hasn’t been explored, or has a FES value that is irrelevanlty high, the cutoff will set the error of that region to zero. If the non_exploration_penalty is larger than zero, the error of that region will take the value of the non_exploration_penalty instead of zero. Default is set to 0.

• use_weighted_st_dev (bool, optional) – When set to True, the calculated error will be the weighted standard deviation ( var^0.5 ). When set to False, the calculated error will be the standard error ( (var/n_sample)^0.5 ). Defaults to True. (The standard devaition is expected to converge after enough time, while the standard error is expected to decrease as more datapoints are added.)

• hp_centre_x (float, optional) – CV1-position of harmonic potential. Defaults to 0.0.

• hp_centre_y (float, optional) – CV2-position of harmonic potential. Defaults to 0.0.

• hp_kappa_x (int, optional) – CV1-force_constant of harmonic potential. Defaults to 0.

• hp_kappa_y (int, optional) – CV2-force_constant of harmonic potential. Defaults to 0.

• lw_centre_x (float, optional) – CV1-position of lower wall potential. Defaults to 0.0.

• lw_centre_y (float, optional) – CV2-position of lower wall potential. Defaults to 0.0.

• lw_kappa_x (int, optional) – CV1-force_constant of lower wall potential. Defaults to 0.

• lw_kappa_y (int, optional) – CV2-force_constant of lower wall potential. Defaults to 0.

• uw_centre_x (float, optional) – CV1-position of upper wall potential. Defaults to 0.0.

• uw_centre_y (float, optional) – CV2-position of upper wall potential. Defaults to 0.0.

• uw_kappa_x (int, optional) – CV1-force_constant of upper wall potential. Defaults to 0.

• uw_kappa_y (int, optional) – CV2-force_constant of upper wall potential. Defaults to 0.

• F_static (array, optional) – Option to provide a starting bias potential that remains constant through the algorithm. This could be a harmonic potential, an previously used MetaD potential or any other bias potential defined on the grid. Defaults to np.zeros((1,1)), which will automatically set F_static to a zero-array with shape=nbins.

• compare_to_reference_FES (int, optional) – Do you wanto to compare with a reference FES (dev option). Defaults to 0.

• ref_fes (array, optional) – Reference FES (dev options). Defaults to an array of zeros.

Returns

[X, Y, Ftot_den, Ftot_x, Ftot_y, ofv, ofe, cutoff, volume_history, ofe_history, time_history, Ftot_den2, ofv_num_x, ofv_num_y]

X (array of size (nbins[1], nbins[0])): CV1 grid positions

Y (array of size (nbins[1], nbins[0])): CV2 grid positions

Ftot_den (array of size (nbins[1], nbins[0])): Cumulative biased probability density, equivalent to an unbiased histogram of samples in CV space.

Ftot_x (array of size (nbins[1], nbins[0])): CV1 component of the Mean Force.

Ftot_y (array of size (nbins[1], nbins[0])): CV2 component of the Mean Force.

ofv (array of size (nbins[1], nbins[0])): on the fly variance estimate of the local convergence

ofe (array of size (nbins[1], nbins[0])): on the fly error estimate of the local convergence

cutoff (array of size (nbins[1], nbins[0])): Array of ones and zeros. By default array of only ones. When the FES or probablity density (Ftot_den) values are outside their respective cutoff, the corresponding cutoff value will be a zero.

volume_history (list of size (nbins[1], nbins[0])): List of scalars indicating the explored volume, as a percentage of the ones in the cutoff array.

ofe_history (list of size (total_number_of_hills/error_pace,)): Running estimate of the global convergence of the mean force by calculating the statistical error of the mean force. Error calculated and added to list with error_pace.

ofe_history_window (optional)(list of size (total_number_of_hills/error_pace)): running estimate of the ofe within the “window” (specified and activated by using the input window_corners).

time_history (list of size (total_number_of_hills/error_pace,)): time array of volume_history, ofe_history and ofe_history_window if applicable.

Ftot_den2 (array of size (nbins[1], nbins[0])): Cumulative squared biased probability density

ofv_x (array of size (nbins[1], nbins[0])): intermediate component in the calculation of the CV1 “on the fly variance” ( sum of: pb_t * dfds_x ** 2)

ofv_y (array of size (nbins[1], nbins[0])): intermediate component in the calculation of the CV2 “on the fly variance” ( sum of: pb_t * dfds_y ** 2)

Return type

list

pyMFI.MFI.bootstrap_2D_new(X, Yrow, force_array, n_bootstrap, min_grid=array([- 3, - 3]), max_grid=array([3, 3]), periodic=array([0, 0]), FES_cutoff=0, Ftot_den_cutoff=0, non_exploration_penalty=0)

Algorithm to determine bootstrap error. Takes in a collection of force-terms and with each itteration, a random selection of force-terms will be used to calculate a FES. The average and st.dev of all FESes will be calculated.

Parameters

• X (array of size (nbins[1], nbins[0])) – CV1 grid positions

• Y (array of size (nbins[1], nbins[0])) – CV2 grid positions

• force_array (list) – collection of force terms (n * [Ftot_den, Ftot_x, Ftot_y])

• n_bootstrap (int) – bootstrap iterations

• min_grid (array, optional) – Lower bound of the force domain. Defaults to np.array((-np.pi, -np.pi)).

• max_grid (array, optional) – Upper bound of the force domain. Defaults to np.array((np.pi, np.pi)).

• periodic (list or array of shape (2,))) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

• FES_cutoff (float, optional) – Cutoff applied to error calculation for FES values over the FES_cutoff. If the cutoff applies, the error will be set to zero, otherwise the error will stay the same. Defaults to 0. When FES_cutoff <= 0, no cufoff is applied. Use with care, computing the fes in the loop renders the calculation slow.

• Ftot_den_cutoff (float, optional) – Cutoff applied to error calculation for probability density (Ftot_den) values below the Ftot_den_cutoff. If the cutoff applies, the error will be set to zero, otherwise the error will stay the same. Defaults to 0.1. When Ftot_den_cutoff <= 0, no cufoff is applied.

• non_exploration_penalty (float, optional) – Turns zero-value error to the non_exploration_penalty value. This should be used in combination with the cutoff. If some part of CV space hasn’t been explored, or has a FES value that is irrelevanlty high, the cutoff will set the error of that region to zero. If the non_exploration_penalty is larger than zero, the error of that region will take the value of the non_exploration_penalty instead of zero. Default is set to 0.

Returns

[FES_avr, sd_fes, sd_fes_prog ]

FES_avr (array of size (nbins[1], nbins[0])): Average of all FESes generated.

sd_fes (array of size (nbins[1], nbins[0])): Map of standard deviation of all FESes generated.

sd_fes_prog (array of size (n_bootstrap,)): The standard deviation of all FESes generated after each itteration.

Return type

list

pyMFI.MFI.bootstrap_2D_to_1D(force_array, n_bootstrap, kT=1, error_pace=- 1, min_grid=array([- 3, - 3]), max_grid=array([3, 3]), periodic=array([0, 0]), Ftot_den_cutoff=0.1, use_weighted_st_dev=True)

Algorithm to determine bootstrap error in x-dimension and the error in y-dimension using force gradients defined in x-y-dimension. Takes in a collection of force-terms and with each itteration, a random selection of force-terms will be used to calculate a FES in both x- and y-dimension. The average and st.dev of all FESes will be calculated.

Parameters

• force_array (list) – collection of force terms (n * [Ftot_den, Ftot_x, Ftot_y])

• n_bootstrap (int) – bootstrap iterations

• kT (float, optional) – Temperature of the system.

• error_pace (int, optional) – Pace at which the error is calculated. Pace of 10 means that the error will be calculated every 10 iterations. Pace of -1 will calculate the error a total of 100 times. Default is set to 0.

• min_grid (array, optional) – Lower bound of the force domain. Defaults to np.array((-3, -3)).

• max_grid (array, optional) – Upper bound of the force domain. Defaults to np.array((3, 3)).

• periodic (list or array of shape (2,))) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

• Ftot_den_cutoff (float, optional) – Cutoff applied to error calculation for probability density (Ftot_den) values below the Ftot_den_cutoff. If the cutoff applies, the error will be set to zero, otherwise the error will stay the same. Defaults to 0.1. When Ftot_den_cutoff <= 0, no cufoff is applied.

• use_weighted_st_dev (bool, optional) – When set to True, the calculated error will be the weighted standard deviation ( var^0.5 ). When set to False, the calculated error will be the standard error ( (var/n_sample)^0.5 ). Defaults to True. (The standard devaition is expected to converge after enough time, while the standard error is expected to decrease as more datapoints are added.)

Returns

[FES_x_avr, error_x, sd_fes_x_prog, FES_y_avr, error_y, sd_fes_y_prog ]

FES_x_avr (array of size (nbins[1], nbins[0])): Average of all FESes generated in x-dimension.

error_x (array of size (nbins[1], nbins[0])): Map of standard deviation of all FESes generated in x-dimension.

sd_fes_x_prog (array of size (n_bootstrap,)): The standard deviation of all FESes generated in x-dimension after each itteration.

FES_y_avr (array of size (nbins[1], nbins[0])): Average of all FESes generated in y-dimension.

error_y (array of size (nbins[1], nbins[0])): Map of standard deviation of all FESes generated in y-dimension.

sd_fes_y_prog (array of size (n_bootstrap,)): The standard deviation of all FESes generated in y-dimension after each itteration.

Return type

list

pyMFI.MFI.coft(HILLS='HILLS', FES='FES', kT=1, WellTempered=- 1, total_number_of_hills=100, stride=10, min_grid=array([- 3.14159265, - 3.14159265]), max_grid=array([3.14159265, 3.14159265]), nbins=array([200, 200]))

Compute a time-independent estimate of the Mean Thermodynamic Force, i.e. the free energy gradient in 2D CV spaces.

Parameters

• HILLS (array) – HILLS array. Defaults to “HILLS”.

• FES (array of size (nbins[0], nbins[1])) –

• kT (float, optional) – Boltzmann constant multiplied with temperature (reduced format, 120K -> 1).

• WellTempered (binary, optional) – Is the simulation well tempered? 1 or yes and 0 for no. Defaults to 1.

• total_number_of_hills (int, optional) – Number of HILLS to analyse. Defaults to 100.

• min_grid (array, optional) – Lower bound of the force domain. Defaults to np.array((-np.pi, -np.pi)).

• max_grid (array, optional) – Upper bound of the force domain. Defaults to np.array((np.pi, np.pi)).

• nbins (array, optional) – number of bins in CV1,CV2. Defaults to np.array((200,200)).

Returns

[c,Bias]

c (array of size (total_numer_of_hills,)): Ensemble average of the bias in the unperturbed ensemble, calculated after the deposition of each metadynamics hill.

Bias (array of size (nbins[0], nbins[1])): Metadynamics Bias reconstructed from HILLS data.

Return type

list

pyMFI.MFI.find_hp_force(hp_centre_x, hp_centre_y, hp_kappa_x, hp_kappa_y, X, Y, min_grid, max_grid, grid_space, periodic)

Find 2D harmonic potential force.

Parameters

• hp_centre_x (float) – CV1-position of harmonic potential

• hp_centre_y (float) – CV2-position of harmonic potential

• hp_kappa_x (float) – CV1-force_constant of harmonic potential

• hp_kappa_y (float) – CV2-force_constant of harmonic potential

• X (array) – 2D array of CV1 grid positions

• Y (array) – 2D array of CV2 grid positions

• min_grid (list) – list of CV1-minimum value of grid and CV2-minimum value of grid

• max_grid (list) – list of CV1-maximum value of grid and CV2-maximum value of grid

• grid_space (list) – list of CV1-grid spacing and CV2-grid spacing

• periodic (list or array of shape (2,))) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

Returns

[F_harmonic_x, F_harmonic_y] F_harmonic_x F_harmonic_y

Return type

list

pyMFI.MFI.find_lw_force(lw_centre_x, lw_centre_y, lw_kappa_x, lw_kappa_y, X, Y, min_grid, max_grid, grid_space, periodic)

Find lower half of 2D harmonic potential force equivalent to f = 2 * lw_kappa * (grid - lw_centre) for grid < lw_centre and f = 0 otherwise. This can change for periodic cases.

Parameters

• lw_centre_x (float) – CV1-position of lower wall potential

• lw_centre_y (float) – CV2-position of lower wall potential

• lw_kappa_x (float) – CV1-force_constant of lower wall potential

• lw_kappa_y (float) – CV2-force_constant of lower wall potential

• X (array) – 2D array of CV1 grid positions

• Y (array) – 2D array of CV2 grid positions

• min_grid (list) – list of CV1-minimum value of grid and CV2-minimum value of grid

• max_grid (list) – list of CV1-maximum value of grid and CV2-maximum value of grid

• grid_space (list) – list of CV1-grid spacing and CV2-grid spacing

• periodic (list or array of shape (2,))) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

Returns

[F_wall_x, F_wall_y] F_wall_x F_wall_y

Return type

list

pyMFI.MFI.find_periodic_point(x_coord, y_coord, min_grid, max_grid, periodic)

Finds periodic copies of input coordinates. First checks if systems is periodic. If not, returns input coordinate array. Next, it checks if each coordinate is within the boundary range (grid min/max +/- grid_ext). If it is, periodic copies will be made on the other side of the CV-domain.

Parameters

• x_coord (float) – CV1-coordinate

• y_coord (float) – CV2-coordinate

• min_grid (list) – list of CV1-minimum value of grid and CV2-minimum value of grid

• max_grid (list) – list of CV1-maximum value of grid and CV2-maximum value of grid

• periodic (list or array of shape (2,)) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

Returns

list of [x-coord, y-coord] pairs (i.e. [[x1,y1], [x2,y2], …, [xn,yn]])

Return type

list

pyMFI.MFI.find_uw_force(uw_centre_x, uw_centre_y, uw_kappa_x, uw_kappa_y, X, Y, min_grid, max_grid, grid_space, periodic)

Find upper half of 2D harmonic potential force equivalent to f = 2 * uw_kappa * (grid - uw_centre) for grid > uw_centre and f = 0 otherwise. This can change for periodic cases.

Parameters

• lw_centre_x (float) – CV1-position of upper wall potential

• lw_centre_y (float) – CV2-position of upper wall potential

• lw_kappa_x (float) – CV1-force_constant of upper wall potential

• lw_kappa_y (float) – CV2-force_constant of upper wall potential

• X (array) – 2D array of CV1 grid positions

• Y (array) – 2D array of CV2 grid positions

• min_grid (list) – list of CV1-minimum value of grid and CV2-minimum value of grid

• max_grid (list) – list of CV1-maximum value of grid and CV2-maximum value of grid

• grid_space (list) – list of CV1-grid spacing and CV2-grid spacing

• periodic (list or array of shape (2,))) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

Returns

[F_wall_x, F_wall_y] F_wall_x F_wall_y

Return type

list

pyMFI.MFI.get_cutoff(Ftot_den=None, Ftot_den_cutoff=0.1, FX=None, FY=None, FES_cutoff=- 1)

Finds the cutoff array according to the specifications.

Parameters

• Ftot_den (np.array, optional) – If a probability density (Ftot_den) cutoff should be applied, this argument in necessary. Defaults to None.

• Ftot_den_cutoff (float, optional) – Specifies the cutoff limit of the probability density. When negative or zero, no probability density cutoff is applied. Defaults to -1.

• FX (np.array, optional) – Force gradient of X or CV1. If a free energy surface (FES) cutoff should be applied, this argument in necessary. Defaults to None.

• FY (np.array, optional) – Force gradient of Y or CV2. If a free energy surface (FES) cutoff should be applied, this argument in necessary. Defaults to None.

• FES_cutoff (list or float, optional) – Required list: [FES_cutoff_limit, min_grid, max_grid, periodic]. If list is not provided, FES_cutoff will not be applied. Defaults to -1.

Returns

cutoff array with the shape of FY. Elements that correspond to the probability density above the Ftot_den_cutoff or the FES below the FES_cutoff will be 1. Elements outside the cutoff will be 0.

Return type

np.array

pyMFI.MFI.get_cutoff_1D(Ftot_den=None, Ftot_den_cutoff=0.1, FES=None, FES_cutoff=23)

Finds the cutoff array according to the specifications.

Parameters

• Ftot_den (np.array, optional) – If a probability density (Ftot_den) cutoff should be applied, this argument in necessary. Defaults to None.

• Ftot_den_cutoff (float, optional) – Specifies the cutoff limit of the probability density. Values below the limit will be cut-off. Will only be applied if Ftot_den is provided. Defaults to 0.1.

• FES (np.array, optional) – free energy surfacee. If a free energy surface (FES) cutoff should be applied, this argument in necessary. Defaults to None.

• FES_cutoff (float, optional) – Specifies the cutoff limit of the FES. Values above the limit will be cut-off. Will only be applied if FES is provided. Defaults to 23.

Returns

cutoff array with the shape of Ftot_den or FES. Elements that correspond to the probability density above the Ftot_den_cutoff or the FES below the FES_cutoff will be 1. Elements outside the cutoff will be 0.

Return type

array

pyMFI.MFI.index(position, min_grid, grid_space)

Finds (approximate) index of a position in a grid. Independent of CV-type.

Parameters

• position (float) – position of interest

• min_grid (float) – minimum value of grid

• grid_space (float) – grid spacing

Returns

index of position

Return type

int

pyMFI.MFI.intg_2D(FX, FY, min_grid=array([- 3.14159265, - 3.14159265]), max_grid=array([3.14159265, 3.14159265]), nbins=array([200, 200]))

2D integration of force gradient (FX, FY) to find FES using finite difference method.

Parameters

• FX (array of size (nbins[1], nbins[0])) – CV1 component of the Mean Force.

• FY (array of size (nbins[1], nbins[0])) – CV2 component of the Mean Force.

• min_grid (array, optional) – Lower bound of the simulation domain. Defaults to np.array((-np.pi, -np.pi)).

• min_grid – Upper bound of the simulation domain. Defaults to np.array((np.pi, np.pi)).

• nbins (int, optional) – number of bins in CV1,CV2. Defaults to np.array((200,200)).

Returns

[X, Y, fes]

X (array of size (nbins[1], nbins[0])): CV1 grid positions

Y (array of size (nbins[1], nbins[0])): CV2 grid positions

fes (array of size (nbins[1], nbins[0])): Free Energy Surface

Return type

list

pyMFI.MFI.intgrad2(fx, fy, min_grid=array([- 2, - 2]), max_grid=array([2, 2]), periodic=array([0, 0]), intconst=0)

This function uses the inverse of the gradient to reconstruct the free energy surface from the mean force components. [John D’Errico (2022). Inverse (integrated) gradient (https://www.mathworks.com/matlabcentral/fileexchange/9734-inverse-integrated-gradient), MATLAB Central File Exchange. Retrieved May 17, 2022.] [Translated from MatLab to Python by Francesco Serse (https://github.com/Fserse)]

Parameters

• fx (array) – (ny by nx) array. X-gradient to be integrated.

• fy (array) – (ny by nx) array. X-gradient to be integrated.

• nx (integer) – nuber of datapoints in x-direction. Default to 0: will copy the shape of the input gradient.

• ny (integer) – nuber of datapoints in y-direction. Default to 0: will copy the shape of the input gradient.

• intconst (float) – Minimum value of output FES

• periodic (list or array of shape (2,))) – Specifies if system is periodic. First entry sepcifies if CV1 is periodic. Second entry sepcifies if CV2 is periodic. value of 0 corresponds to non-periodic CV1. Value of 1 corresponds to periodic CV1.

• min_grid (list/array of length=2) – list/array of minimum value of [x-grid, y-grid]

• max_grid (list/array of length=2) – list/array of maximum value of [x-grid, y-grid]

• nbins (list/array of length=2) – list/array of number of data pointis of [x-grid, y-grid]. Default to 0: will copy the shape of the input gradient.

Returns

[X, Y, fhat]

X (ny by nx array): X-component of meshgrid

Y (ny by nx array): Y-component of meshgrid

fhat (ny by nx array): integrated free energy surface

Return type

list

pyMFI.MFI.load_HILLS_2D(hills_name='HILLS')

Load 2-dimensional hills data (includes time, position_x, position_y, hills_parameters ).

Parameters

hills_name (str, optional) – Name of hills file. Defaults to “HILLS”.

Returns

hills data with length equal to the total number of hills. Information: [time [ps], position_x [nm], position_y [nm], MetaD_sigma_x [nm], MetaD_sigma_y [nm], MetaD_height [nm], MetaD_biasfactor]

Return type

np.array

pyMFI.MFI.load_npy(name)

Loads np.array of a file with .npy format

Parameters

name (string) – Name of file

Returns

object to be loaded. Must be a numpy array.

Return type

np.array

pyMFI.MFI.load_pkl(name)

Loads list/array of a file with .pkl format

Parameters

name (string) – Name of file

Returns

object to be loaded

Return type

any

pyMFI.MFI.load_position_2D(position_name='position')

Load 2-dimensional position/trajectory data.

Parameters

position_name (str, optional) – Name of position file. Defaults to “position”.

Returns

[position_x, position_y] position_x (np.array of shape (number_of_positions,)): position (or COLVAR) data of x-dimension (or CV1) position_y (np.array of shape (number_of_positions,)): position (or COLVAR) data of y-dimension (or CV2)

Return type

list

pyMFI.MFI.mean_force_variance(Ftot_den, Ftot_den2, Ftot_x, Ftot_y, ofv_num_x, ofv_num_y, Ftot_den_limit=1e-10, use_weighted_st_dev=True)

Calculates the variance of the mean force

Parameters

• Ftot_den (array of size (nbins[1], nbins[0])) – Cumulative biased probability density

• 25Ftot_den2 (array of size (nbins[1], nbins[0])) – Cumulative squared biased probability density

• Ftot_x (array of size (nbins[1], nbins[0])) – CV1 component of the Mean Force.

• Ftot_y (array of size (nbins[1], nbins[0])) – CV2 component of the Mean Force.

• ofv_num_x (array of size (nbins[1], nbins[0])) – intermediate component in the calculation of the CV1 “on the fly variance” ( sum of: pb_t * dfds_x ** 2)

• ofv_num_y (array of size (nbins[1], nbins[0])) – intermediate component in the calculation of the CV2 “on the fly variance” ( sum of: pb_t * dfds_y ** 2)

• Ftot_den_limit (scalar) – threshold in the cumulative biased probability density where data is discarded. Defaults to 0.

Returns

[ofv, ofe]

ofv (array of size (nbins[1], nbins[0])): “on the fly variance”

ofe (array of size (nbins[1], nbins[0])): “on the fly error”

Return type

list

pyMFI.MFI.patch_2D(master_array)

Takes in a collection of force terms and patches them togehter to return the patched force terms

Parameters

master_array (list) – collection of force terms (n * [Ftot_den, Ftot_den2, Ftot_x, Ftot_y, ofv_x, ofv_y])

Returns

[FP, FP2, FX, FY, OFV_X, OFV_Y]

FP(array of size (nbins[1], nbins[0])): Patched probability density

FP2(array of size (nbins[1], nbins[0])): Patched (probability density squared)

FX(array of size (nbins[1], nbins[0])): Patched Froce gradient in x-direction (CV1 direction)

FY(array of size (nbins[1], nbins[0])): Patched Froce gradient in y-direction (CV2 direction)

OFV_X (array of size (nbins[1], nbins[0])): “on the fly variance”-term for the calculation of the variance of the froce gradient in x-direction (CV1 direction)

OFV_Y (array of size (nbins[1], nbins[0])): “on the fly variance”-term for the calculation of the variance of the froce gradient in y-direction (CV2 direction)

Return type

list

pyMFI.MFI.patch_2D_simple(master_array)

Takes in a collection of force and patches only the probability density and mean forces

Parameters

master_array (list) – collection of force terms (n * [Ftot_den, Ftot_x, Ftot_y])

Returns

[FP, FP2, FX, FY, OFV_X, OFV_Y]

FP (array of size (nbins[1], nbins[0])): Patched probability density

FX (array of size (nbins[1], nbins[0])): Patched Froce gradient in x-direction (CV1 direction)

FY (array of size (nbins[1], nbins[0])): Patched Froce gradient in y-direction (CV2 direction)

Return type

list

pyMFI.MFI.patch_to_base_variance(master0, master, Ftot_den_limit=1e-10, use_weighted_st_dev=True)

Patches force terms of a base simulation (alaysed prior to current simulation) with current simulation to return patched “on the fly variance”.

Parameters

• master0 (list) – Force terms of base simulation (alaysed prior to current simulation) [Ftot_den, Ftot_den2, Ftot_x, Ftot_y, ofv_num_x, ofv_num_y]

• master (list) – Force terms of current simulation [Ftot_den, Ftot_den2, Ftot_x, Ftot_y, ofv_num_x, ofv_num_y]

• Ftot_den_limit (int) – Truncates the probability density below Ftot_den_limit. Default set to 10**-10.

• use_weighted_st_dev (bool) –

Returns

[PD_patch, FX_patch, FY_patch, OFV, OFE]

PD_patch (array of size (nbins[1], nbins[0])): Patched probability density

FX_patch (array of size (nbins[1], nbins[0])): Patched Froce gradient in x-direction (CV1 direction)

FY_patch (array of size (nbins[1], nbins[0])): Patched Froce gradient in y-direction (CV2 direction)

OFV (array of size (nbins[1], nbins[0])): modulus of patched “on the fly variance”

OFE(array of size (nbins[1], nbins[0])): modulus of patched “on the fly error”

Return type

list

pyMFI.MFI.plot_bootstrap(X, Y, FES, sd_fes, sd_fes_prog, FES_lim=11, sd_lim=11, FES_levels=None, sd_levels=None)

Plots result of bootstrap analysis. 1. Average FES, 2. average varinace, 3. variance progression

Parameters

• X (array of size (nbins[1], nbins[0])) – CV1 grid positions

• Y (array of size (nbins[1], nbins[0])) – CV2 grid positions

• FES (array of size (nbins[1], nbins[0])) – Free Energy Surface

• sd_fes (array of size (nbins[1], nbins[0])) – Map of standard deviation of all FESes generated.

• sd_fes_prog (list / np.array of size (bootstrap_iterations,)) – Progression of the standard deviation (of all FESes generated after each bootstrap itteration).

• FES_lim (int, optional) – Upper energy limit of FES plot. Defaults to 11.

• sd_lim (int, optional) – Upper variance limit of variance plot. Defaults to 11.

• FES_levels (int, optional) – Amout of contour levels shown in FES plot. Default is set to None, in which case FES_levels = int(FES_lim + 1).

• ofe_levels (int, optional) – Amout of contour levels shown in standard deviation plot. Default is set to None, in which case FES_levels = int(FES_lim + 1).

pyMFI.MFI.plot_bootstrap_2D_to_1D(X, Y, FES_2D, FES_x, FES_y, error_x, error_y, FES_lim=20)

pyMFI.MFI.plot_patch_2D(X, Y, FES, TOTAL_DENSITY, lim=50)

Plots 1. FES, 2. Cumulative biased probability density

Parameters

• X (array of size (nbins[1], nbins[0])) – CV1 grid positions

• Y (array of size (nbins[1], nbins[0])) – CV2 grid positions

• FES (array of size (nbins[1], nbins[0])) – Free Energy Surface

• TOTAL_DENSITY (array of size (nbins[1], nbins[0])) – Cumulative biased probability density

pyMFI.MFI.plot_recap_2D(X, Y, FES, TOTAL_DENSITY, CONVMAP, CONV_history, CONV_history_time, FES_lim=50, ofe_map_lim=50, FES_step=1, ofe_step=1)

Plots 1. FES, 2. varinace_map, 3. Cumulative biased probability density, 4. Convergece of variance.

Parameters

• X (array of size (nbins[1], nbins[0])) – CV1 grid positions

• Y (array of size (nbins[1], nbins[0])) – CV2 grid positions

• FES (array of size (nbins[1], nbins[0])) – Free Energy Surface

• TOTAL_DENSITY (array of size (nbins[1], nbins[0])) – Cumulative biased probability density

• CONVMAP (array of size (nbins[1], nbins[0])) – varinace_map

• CONV_history (list) – Convergece of variance

• CONV_history_time (list) – Simulation time corresponding to CONV_history

pyMFI.MFI.print_progress(iteration, total, bar_length=50, variable_name='progress variable', variable=0)

Function to show a progress bar, that fills up as the iteration number reaches the total. Prints a variable at the end of the progress bar, that can be continiously updated.

Parameters

• iteration (int) – Currrent iteration

• total (int) – Total iterations

• bar_length (int, optional) – Length of the progress bar. Defaults to 50.

• variable_name (str, optional) – Name of variable that is being shown at the end of the progress bar. Defaults to ‘progress variable’.

• variable (float, optional) – Variable that is being shown at the end of the progress bar. Defaults to 0.

pyMFI.MFI.reduce_to_window(input_array, min_grid, grid_space, x_min=- 0.5, x_max=0.5, y_min=- 1.5, y_max=1.5)

Reduces an 2D input array to a specified range.

Parameters

• input_array (array) – 2D array to be reduced

• min_grid (list) – list of CV1-minimum value of grid and CV2-minimum value of grid

• grid_space (list) – list of CV1-grid spacing and CV2-grid spacing

• x_min (float, optional) – lower CV1-value of output array. Defaults to -0.5.

• x_max (float, optional) – upper CV1-value of output array. Defaults to 0.5.

• y_min (float, optional) – lower CV2-value of output array. Defaults to -1.5.

• y_max (float, optional) – upper CV2-value of output array. Defaults to 1.5.

Returns

reduced array

Return type

array

pyMFI.MFI.save_npy(object, file_name)

Saves np.array in a file with .npy format

Parameters

• object (np.array) – object to be saved. Must be a numpy array.

• file_name (string) – Name of file

pyMFI.MFI.save_pkl(object, file_name)

Saves a list/array in a file with .pkl format

Parameters

• object (any) – object to be saved

• file_name (string) – Name of file

pyMFI.MFI.zero_to_nan(input_array)

Function to turn all zero-elements to np.nan. Works for any shapes.

Parameters

input_array (array of arbitrary shape) – non specific array

Returns

input array with zero-elements turned to np.nan

Return type

array

pyMFI.MFI1D module

pyMFI.MFI1D.MFI_1D(HILLS, position, bw=0.1, kT=1, min_grid=- 2, max_grid=2, nbins=201, log_pace=- 1, error_pace=- 1, WellTempered=1, nhills=- 1, periodic=0, hp_centre=0.0, hp_kappa=0, lw_centre=0.0, lw_kappa=0, uw_centre=0.0, uw_kappa=0, F_static=array([0.]), Ftot_den_limit=1e-10, FES_cutoff=0, Ftot_den_cutoff=0, non_exploration_penalty=0, save_intermediate_fes_error_cutoff=False, use_weighted_st_dev=True)

Compute a time-independent estimate of the Mean Thermodynamic Force, i.e. the free energy gradient in 1D CV spaces.

Parameters

• HILLS (array) – HILLS array from HILLS file from metadynamics simulation using PLUMED. Contains number_of_hills * [time, position, metaD_width, metaD_height metaD_biasfactor].

• position (str) – CV/position array. Defaults to “position”.

• bw (float, optional) – bandwidth for the construction of the KDE estimate of the biased probability density. Defaults to 1.

• kT (float, optional) – Boltzmann constant multiplied with temperature (reduced format, 120K -> 1).

• min_grid (float, optional) – Lower bound of the force domain. Defaults to -2.

• max_grid (float, optional) – Upper bound of the force domain. Defaults to 2.

• nbins (int, optional) – number of bins in grid. Defaults to 101.

• log_pace (int, optional) – Pace for outputting progress and convergence. Defaults to -1. When set to -1, progress will be outputted 5 times in total.

• error_pace (int, optional) – Pace for calculating the on-the-fly error for estimating global convergence. Defaults to -1. When set to -1, on-the-fly error will be calculated 50 times.

• WellTempered (binary, optional) – Is the simulation well tempered?. Defaults to 1.

• nhills (binary, optional) – Number of HILLS to be analysed. When set to -1, all HILLS will be analysed. Defaults to -1.

• periodic (binary, optional) – Is the CV space periodic? 1 for yes. Defaults to 0.

• hp_centre (float, optional) – position of harmonic potential. Defaults to 0.0.

• hp_kappa (flaot, optional) – force_constant of harmonic potential. Defaults to 0.

• lw_centre (float, optional) – position of lower wall potential. Defaults to 0.0.

• lw_kappa (flaot, optional) – force_constant of lower wall potential. Defaults to 0.

• uw_centre (float, optional) – position of upper wall potential. Defaults to 0.0.

• uw_kappa (flaot, optional) – force_constant of upper wall potential. Defaults to 0.

• F_static (array, optional) – Option to provide a starting bias potential that remains constant through the algorithm. This could be a harmonic potential, an previously used MetaD potential or any other bias potential defined on the grid. Defaults to np.zeros(1), which will automatically set F_static to a zero-array with length=nbins.

• Ftot_den_limit (float, optional) – Probability density limit below which data will be set to zero (this is done for numerical stability reasons. For default Ftot_den_limit, numerical difference is negligable). Defaults to 1E-10.

• FES_cutoff (float, optional) – Cutoff applied to FES and error calculation for FES values over the FES_cutoff. All FES values above the FES_cutoff won’t contribuite towards the error. Useful when high FES values have no physical meaning. When FES_cutoff = 0, no cufoff is applied.Defaults to 0.

• Ftot_den_cutoff (flaot, optional) – Cutoff applied to FES and error calculation for Ftot_den (Probability density) values below the Ftot_den_cutoff. All FES values that are excluded by the cutoff won’t contribuite towards the error. Useful when low Probability density values have little statistical significance or no physical meaning. When Ftot_den_cutoff = 0, no cufoff is applied. Defaults to 0. non_exploration_penalty (float, optional): Turns zero-value error to the non_exploration_penalty value. This should be used in combination with the cutoff. If some part of CV space hasn’t been explored, or has a FES value that is irrelevanlty high, the cutoff will set the error of that region to zero. If the non_exploration_penalty is larger than zero, the error of that region will take the value of the non_exploration_penalty instead of zero. Default is set to 0.

• save_intermediate_fes_error_cutoff (bool, optional) – If Ture, every time the error is calculated, the FES, variance, standard deviation and cutoff will be saved. Defaults to False.

• use_weighted_st_dev (bool, optional) – When set to True, the calculated error will be the weighted standard deviation ( var^0.5 ). When set to False, the calculated error will be the standard error ( (var/n_sample)^0.5 ). Defaults to True. (The standard devaition is expected to converge after enough time, while the standard error is expected to decrease as more datapoints are added.)

Returns

grid, Ftot_den, Ftot_den2, Ftot, ofv_num, FES, ofv, ofe, cutoff, error_evol, fes_error_cutoff_evol

grid (array of size (nbins,)): CV-array.

Ftot_den (array of size (nbins,)): Cumulative biased probability density.

Ftot_den2 (array of size (nbins,): Cumulative (biased probability density squared). Used for error calculation.

Ftot (array of size (nbins,)): Local Mean Force. When integrated will be the FES.

ofv_num (array of size (nbins,)): Numerator of “on-the-fly” variance (ofv). Used for error calculation.

FES (array of size (nbins,)): Free Energy Surface

ofv (array of size (nbins,)): “on-the-fly” variance

ofe (array of size (nbins,)): “on-the-fly” estimate of the standard deviation of the mean force

cutoff (binary array of size (nbins,)): If FES_cutoff and/or Ftot_den_cutoff are specified, grid values outside the cufoff will be zero and grid values inside the cutoff will be one. If FES_cutoff and Ftot_den_cutoff not active, cutoff will be an array of ones.

error_evol (array of size (4, total_number_of_hills//error_pace)): Evolution of error. First array is the collection of the global (average) “on-the-fly” variance, second array is the collection of the global (average) “on-the-fly” standard deviation, third array is the percentage of values within the specified cutoff, and the fourth array are the simulation times of the latter former arrays.

fes_error_cutoff_evol (array of size (4, total_number_of_hills//error_pace, nbins)): First array is the collection of FES, second array is the collection of local “on-the-fly” variance, third array is the collection of local “on-the-fly” standard deviation, and the fourth is the collection of the cutoffs. The elements of the collections were calculated at the simulation times of the fourth array in error_evol.

Return type

tuple

pyMFI.MFI1D.bootstrap_1D(grid, force_array, n_bootstrap, set_fes_minima=None)

Algorithm to determine bootstrap error

Parameters

• grid (array of shape (nbins,)) – CV grid positions

• force_array (list) – collection of force terms (n * [Ftot_den, Ftot]).

• n_bootstrap (int) – bootstrap iterations

Returns

Average of all FES generated during the bootstrap algorithm. sd_fes (array of shape (nbins,)): Standard deviation of all FES generated during the bootstrap algorithm. sd_fes_prog (array of shape (n_bootstrap,)): Global average of the standard deviation after each bootstrap iteration. When this array converges, enough iterations have been performed. If it does not converge, move iterations are necessary.

Return type

FES_avr (array of shape (nbins,))

pyMFI.MFI1D.bootstrap_forw_back(grid, forward_force, backward_force, n_bootstrap, set_fes_minima=None)

Algorithm to determine bootstrap error

Parameters

• grid (array of shape (nbins,)) – CV grid positions

• forward_force (list) – collection of force terms (n * [Ftot_den, Ftot]) from forward transition

• backward_force (list) – collection of force terms (n * [Ftot_den, Ftot]) from backward transition

• n_bootstrap (int) – bootstrap iterations

• set_fes_minima (str, optional) – USed to specify how to set the minima of the FES. When set to “first”, the first element of the fes will be set to 0, and the rest of the FES array will be shifted by the same amount. When set to None, the smalles element of the FES will be set to 0, and the rest of the FES array will be shifted by the same amount. Defauts is None.

Returns

[FES_avr, sd_fes, sd_fes_prog]

FES_avr (array of shape (nbins,)): Average of all FES generated during the bootstrap algorithm.

sd_fes (array of shape (nbins,)): Standard deviation of all FES generated during the bootstrap algorithm.

sd_fes_prog (array of shape (n_bootstrap,)): Global average of the standard deviation after each bootstrap iteration. When this array converges, enough iterations have been performed. If it does not converge, move iterations are necessary.

Return type

list

pyMFI.MFI1D.find_hp_force(hp_centre, hp_kappa, grid, min_grid, max_grid, grid_space, periodic)

Find 1D harmonic potential force equivalent to f = hp_kappa * (grid - hp_centre).

Parameters

• hp_centre (float) – position of harmonic potential

• hp_kappa (float) – force_constant of harmonic potential

• grid (array) – CV grid positions

• min_grid (float) – minimum value of grid

• max_grid (float) – maximum value of grid

• grid_space (float) – space between two consecutive grid values

• periodic (binary) – information if system is periodic. value of 0 corresponds to non-periodic system. Value of 1 corresponds to periodic system.

Returns

harmonic force array

Return type

array

pyMFI.MFI1D.find_lw_force(lw_centre, lw_kappa, grid, min_grid, max_grid, grid_space, periodic)

Find lower half of 1D harmonic potential force equivalent to f = 2 * lw_kappa * (grid - lw_centre) for grid < lw_centre and f = 0 otherwise. This can change for periodic cases.

Parameters

• lw_centre (float) – position of lower wall potential

• lw_kappa (float) – force_constant of lower wall potential

• grid (array) – CV grid positions

• min_grid (float) – minimum value of grid

• max_grid (float) – maximum value of grid

• grid_space (float) – space between two consecutive grid values

• periodic (binary) – information if system is periodic. value of 0 corresponds to non-periodic system. Value of 1 corresponds to periodic system.

Returns

lower wall force array

Return type

array

pyMFI.MFI1D.find_periodic_point(x_coord, min_grid, max_grid, periodic, grid_ext)

Finds periodic copies of input coordinates. First checks if systems is periodic. If not, returns input coordinate array. Next, it checks if each coordinate is within the boundary range (grid min/max +/- grid_ext). If it is, periodic copies will be made on the other side of the CV-domain.

Parameters

• x_coord (float) – CV-coordinate

• min_grid (float) – minimum value of grid

• max_grid (float) – maximum value of grid

• periodic (binary) – information if system is periodic. value of 0 corresponds to non-periodic system; function will only return input coordinates. Value of 1 corresponds to periodic system; function will return input coordinates with periodic copies.

• grid_ext (float) – how far outside the domain periodic copies are searched for. E.g. if the domain is from -2 to +2 and the grid ext(end) is set to 1, then periodic copies are only found for input_coordiante < -2 + 1 or input_coordiante > +2 - 1.

Returns

list of input coord and if applicable periodic copies

Return type

list

pyMFI.MFI1D.find_periodic_point_numpy(coord_array, min_grid, max_grid, periodic, grid_ext)

Finds periodic copies of input coordinates. First checks if systems is periodic. If not, returns input coordinate array. Next, it checks if each coordinate is within the boundary range (grid min/max +/- grid_ext). If it is, periodic copies will be made on the other side of the CV-domain.

Parameters

• coord_array (array) – array of CV-coordinate s

• min_grid (float) – minimum value of grid

• max_grid (float) – maximum value of grid

• periodic (binary) – information if system is periodic. value of 0 corresponds to non-periodic system; function will only return input coordinates. Value of 1 corresponds to periodic system; function will return input coordinates with periodic copies.

• grid_ext (float) – how far outside the domain periodic copies are searched for. E.g. if the domain is from -2 to +2 and the grid ext(end) is set to 1, then periodic copies are only found for input_coordiante < -2 + 1 or input_coordiante > +2 - 1.

Returns

list of input coord and if applicable periodic copies

Return type

np.array

pyMFI.MFI1D.find_uw_force(uw_centre, uw_kappa, grid, min_grid, max_grid, grid_space, periodic)

Find upper half of 1D harmonic potential force equivalent to f = 2 * uw_kappa * (grid - uw_centre) for grid > uw_centre and f = 0 otherwise. This can change for periodic cases.

Parameters

• uw_centre (float) – position of upper wall potential

• uw_kappa (float) – force_constant of upper wall potential

• grid (_type_) – CV grid positions

• min_grid (float) – minimum value of grid

• max_grid (float) – maximum value of grid

• grid_space (float) – space between two consecutive grid values

• periodic (binary) – information if system is periodic. value of 0 corresponds to non-periodic system. Value of 1 corresponds to periodic system.

Returns

upper wall force array

Return type

array

pyMFI.MFI1D.get_cutoff_1D(Ftot_den=None, Ftot_den_cutoff=0.1, FES=None, FES_cutoff=- 1)

Finds the cutoff array according to the specifications. If the “@njit” is deactivated, use get_cutoff_1D_nonjit() instead.

Parameters

• Ftot_den (np.array, optional) – If a probability density (Ftot_den) cutoff should be applied, this argument in necessary. Defaults to None.

• Ftot_den_cutoff (float, optional) – Specifies the cutoff limit of the probability density. Values below the limit will be cut-off. Will only be applied if Ftot_den is provided. Defaults to 0.1.

• FES (np.array, optional) – free energy surfacee. If a free energy surface (FES) cutoff should be applied, this argument in necessary. Defaults to None.

• FES_cutoff (float, optional) – Specifies the cutoff limit of the FES. Values above the limit will be cut-off. Will only be applied if FES is provided. Defaults to 23.

Returns

cutoff array with the shape of Ftot_den or FES. Elements that correspond to the probability density above the Ftot_den_cutoff or the FES below the FES_cutoff will be 1. Elements outside the cutoff will be 0.

Return type

array

pyMFI.MFI1D.get_cutoff_1D_nonjit(Ftot_den=None, Ftot_den_cutoff=0.1, FES=None, FES_cutoff=23)

Finds the cutoff array according to the specifications.

Parameters

• Ftot_den (np.array, optional) – If a probability density (Ftot_den) cutoff should be applied, this argument in necessary. Defaults to None.

• Ftot_den_cutoff (float, optional) – Specifies the cutoff limit of the probability density. Values below the limit will be cut-off. Will only be applied if Ftot_den is provided. Defaults to 0.1.

• FES (np.array, optional) – free energy surfacee. If a free energy surface (FES) cutoff should be applied, this argument in necessary. Defaults to None.

• FES_cutoff (float, optional) – Specifies the cutoff limit of the FES. Values above the limit will be cut-off. Will only be applied if FES is provided. Defaults to 23.

Returns

cutoff array with the shape of Ftot_den or FES. Elements that correspond to the probability density above the Ftot_den_cutoff or the FES below the FES_cutoff will be 1. Elements outside the cutoff will be 0.

Return type

array

pyMFI.MFI1D.index(position, min_grid, grid_space)

Finds (approximate) index of a position in a grid. Independent of CV-type.

Parameters

• position (float) – position of interest

• min_grid (float) – minimum value of grid

• grid_space (float) – grid spacing

Returns

index of position

Return type

int

pyMFI.MFI1D.intg_1D(Force, dx)

Integration of 1D gradient using finite difference method (simpson’s method).

Parameters

• Force (array) – Mean force

• dx (float) – grid spacing (i.e. space between consecutive grid entries)

Returns

Free energy surface

Return type

array

pyMFI.MFI1D.load_HILLS(hills_name='HILLS')

Load 1-dimensional hills data (includes time, position_x, position_y, hills_parameters)

Adjustment 2 - Last hill is removed from output array, as no complete set of COLVAR data is available of that hill.

Parameters

hills_name (str, optional) – Name of hills file. Defaults to “HILLS”.

Returns

Hills data with length equal to the total number of hills. Information: [time [ps], position [nm], MetaD_sigma [nm], MetaD_height [nm], MetaD_biasfactor]

Return type

np.array

pyMFI.MFI1D.load_npy(name)

Loads np.array of a file with .npy format

Parameters

name (string) – Name of file

Returns

object to be loaded. Must be a numpy array.

Return type

np.array

pyMFI.MFI1D.load_pkl(name)

Loads list/array of a file with .pkl format

Parameters

name (string) – Name of file

Returns

object to be loaded

Return type

any

pyMFI.MFI1D.load_position(position_name='position')

Load 1-dimensional position (CV) data.

Parameters

position_name (str, optional) – Name of position file. Defaults to “position”.

Returns

position data

Return type

np.array

pyMFI.MFI1D.patch_forces(force_vector, PD_limit=1e-10)

Takes in a collection of force and probability density and patches them. :param force_vector: collection of force terms (n * [Ftot_den, Ftot]) :type force_vector: list :param PD_limit: Probability density limit below which data will be set to zero (this is done for numerical stability reasons. For default Ftot_den_limit, numerical difference is negligable). Defaults to 1E-10. :type PD_limit: float, optional

Returns

[PD_patch, F_patch]

PD_patch (array): Patched probability density

F_patch (array): Patched mean force

Return type

list

pyMFI.MFI1D.patch_forces_ofe(force_vector, PD_limit=1e-10, use_weighted_st_dev=True, ofe_progression=True, Ftot_den_cutoff=- 1, non_exploration_penalty=0)

Takes in an array of force terms, patches them together and calculates the error. The force vector has to be in the format [Ftot_den, Ftot_den2, Ftot, ofv_num].

Parameters

• force_vector (array of shape (n,4, nbins)) – The force terms are tipically outputted from the MFI_1D functions. They should be in the format: [Ftot_den (Probability density), Ftot_den2, Ftot (Mean Force), ofv_num (numerator of on the fly variance)].

• PD_limit (flaot, optional) – Probability density values below the PD_limit will be approximated to be 0, to ensure the numerical stability of the algorithm. Defaults to 1E-10.

• ofe_non_exploration_penalty (flaot, optional) – On-the-fly error values that that are ouside of the PD_cutoff will have a error equal to the ofe_non_exploration_penalty. Absolute deviation values that that are ouside of the PD_cutoff will have a error equal to the ofe_non_exploration_penalty/10. Defaults to 100.

• use_weighted_st_dev (bool, optional) – When set to True, the calculated error will be the weighted standard deviation ( var^0.5 ). When set to False, the calculated error will be the standard error ( (var/n_sample)^0.5 ). Defaults to True. (The standard devaition is expected to converge after enough time, while the standard error is expected to decrease as more datapoints are added.).

• ofe_progression (bool, optional) – If True, the error progression will be returned. If false, error progression will be an array will all elements being the final error. Default set to False. Ftot_den_cutoff (float, optional): Specifies the cutoff limit of the probability density. Values below the limit will be cut-off. Will only be applied if Ftot_den_cutoff is larger than 0. Defaults to -1. non_exploration_penalty (float, optional): Turns zero-value error to the non_exploration_penalty value. This should be used in combination with the cutoff. If some part of CV space hasn’t been explored, or has a FES value that is irrelevanlty high, the cutoff will set the error of that region to zero. If the non_exploration_penalty is larger than zero, the error of that region will take the value of the non_exploration_penalty instead of zero. Default is set to 0.

Returns

PD_patch, PD2_patch, F_patch, OFV_num_patch, OFE, Aofe_progression

PD_patch (array of size (nbins,)): Patched Probability density

PD2_patch (array of size (nbins,)): Patched (Probability density squared)

F_patch (array of size (nbins,)): Patched Mean Force

OFV_num_patch (array of size (nbins,)): Patched numerator of “on the fly” variance

OFE (array of size (nbins,)): On-the-fly error map. Statistical error of the mean force

Aofe_progression (float): List of average OFE. If ofe_progression=False, all values will be the last calcuated average OFE. If ofe_progression=True, the values correspond to the average OFE after each patch.

Return type

tuple

pyMFI.MFI1D.patch_forces_ofe_AD(force_vector, grid, y, PD_limit=1e-10, use_weighted_st_dev=True, error_progression=True, Ftot_den_cutoff=- 1, FES_cutoff=- 1, non_exploration_penalty=0, set_fes_minima=None)

Takes in an array of force terms, patches them together and calculates the error, FES and average deviation. The force vector has to be in the format [Ftot_den, Ftot_den2, Ftot, ofv_num].

Parameters

• force_vector (array of shape (n,4, nbins)) – The force terms are tipically outputted from the MFI_1D functions. They should be in the format: [Ftot_den (Probability density), Ftot_den2, Ftot (Mean Force), ofv_num (numerator of on the fly variance)].

• grid (array of shape (nbins, )) – Grid that CV-space is defined on.

• y (array of shape (nbins, )) – Reference surface that is used to calculate the absolute devaition.

• PD_limit (flaot, optional) – Probability density values below the PD_limit will be approximated to be 0, to ensure the numerical stability of the algorithm. Defaults to 1E-10.

• use_weighted_st_dev (bool, optional) – When set to True, the calculated error will be the weighted standard deviation ( var^0.5 ). When set to False, the calculated error will be the standard error ( (var/n_sample)^0.5 ). Defaults to True. (The standard devaition is expected to converge after enough time, while the standard error is expected to decrease as more datapoints are added.).

• error_progression (bool, optional) – When set to True, the error (AOFE, AAD, volume ratio) will be calculated after every patch. When set to False, all forces are patched and the error is calculated only at the end. Default set to True. Ftot_den_cutoff (float, optional): Specifies the cutoff limit of the probability density. Values below the limit will be cut-off. Will only be applied if Ftot_den_cutoff is larger than 0. Defaults to -1. FES_cutoff (float, optional): Specifies the cutoff limit of the FES. Values above the limit will be cut-off. Will only be applied if FES is provided. Defaults to 23. non_exploration_penalty (float, optional): Turns zero-value error to the non_exploration_penalty value. This should be used in combination with the cutoff. If some part of CV space hasn’t been explored, or has a FES value that is irrelevanlty high, the cutoff will set the error of that region to zero. If the non_exploration_penalty is larger than zero, the error of that region will take the value of the non_exploration_penalty instead of zero. Default is set to 0.

• set_fes_minima (str, optional) – USed to specify how to set the minima of the FES. When set to “first”, the first element of the fes will be set to 0, and the rest of the FES array will be shifted by the same amount. When set to None, the smalles element of the FES will be set to 0, and the rest of the FES array will be shifted by the same amount. Defauts is None.

Returns

PD_patch, PD2_patch, F_patch, OFV_num_patch, FES, AD, aad_progression, OFE, Aofe_progression, volume_progression

PD_patch (array of size (nbins,)): Patched Probability density

PD2_patch (array of size (nbins,)): Patched (Probability density squared)

F_patch (array of size (nbins,)): Patched Mean Force

OFV_num_patch (array of size (nbins,)): Patched numerator of “on the fly” variance

FES (array of size (nbins,)): Free Energy Surface from patched Mean Force

AD (array of size (nbins,)): Absulute Deviation from FES to reference surface

aad_progression (list): List of average AAD. If error_progression=False, all values will be the last calcuated average AAD. If error_progression=True, the values correspond to the average AAD after each patch.

OFE (array of size (nbins,)): On-the-fly error. Statistical error of the mean force

Aofe_progression (list): List of average OFE. If error_progression=False, all values will be the last calcuated average OFE. If error_progression=True, the values correspond to the average OFE after each patch.

volume_progression (list): List of average volume ratio (Percentage of the grid that is not cutoff, e.g. with PD>0.1). If error_progression=False, all values will be the last calcuated average volume ratio. If error_progression=True, the values correspond to the average volume ratio after each patch.

Return type

tuple

pyMFI.MFI1D.plot_recap(X, FES, Ftot_den, ofe, ofe_history, time_history, y_ref=None, FES_lim=40, ofe_lim=10, error_log_scale=1, cv='CV', time_axis_name='simulation time')

Plot result of 1D MFI algorithm. 1. FES, 2. varinace_map, 3. Cumulative biased probability density, 4. Convergece of variance.

Parameters

• X (array) – gird

• FES (array) – Free energy surface

• Ftot_den (array) – _description_

• ofe (array) – Cumulative biased probability density

• ofe_history (list) – on the fly estimate of the local convergence

• time_history (list) – List of time or simulation number corresponding to error history.

• FES_lim (int, optional) – Upper energy value in FES plot. Defaults to 40.

• ofe_lim (int, optional) – Upper error value in FES plot. Defaults to 10.

• error_log_scale (boolean, optional) – Option to make error_conversion plot with a log scale. 1 for log scale. Defaults to 1.

pyMFI.MFI1D.print_progress(iteration, total, bar_length=50, variable_name='progress variable', variable=0)

pyMFI.MFI1D.save_npy(object, file_name)

Saves np.array in a file with .npy format

Parameters

• object (np.array) – object to be saved. Must be a numpy array.

• file_name (string) – Name of file

pyMFI.MFI1D.save_pkl(object, file_name)

Saves a list/array in a file with .pkl format

Parameters

• object (any) – object to be saved

• file_name (string) – Name of file

pyMFI.MFI1D.window_forces(periodic_positions, periodic_hills, grid, sigma_meta2, height_meta, kT, const, bw2, Ftot_den_limit=1e-10)

Takes in two arrays of positions. The periodic_positions are collected from the COLVAR file during a period of constant bias and calculates the force component associated with the probability density. periodic_hills are also positions collected from the HILLS file and calculates the force component resulting from the metadynamics bias. In a periodic system, positions have periodic copies if applicable.

Parameters

• periodic_positions (array of shape (n,)) – This array contains the positions from the COLVAR file, that were collected during a period of constant bias. n is the number of positions including their periodic copies if applicable.

• periodic_hills (array of shape (n,)) – This array contains the position of the metadynamics hill (from the HILLS file), that was deposited at the beginning of the period of constant bias. If applicable the array also contains periodic copies.

• grid (array of shape (nbins,)) – CV-array.

• sigma_meta2 (float) – width of metadynamics hill squared.

• height_meta (float) – height of metadynamics hill

• kT (float) – Temperature of the system.

• const (float) – constant factor used for the calculation of the kernel. This constant enables the patching of data with different values of bw and stride . const = (1 / (bw * np.sqrt(2 * np.pi) * stride))

• bw2 (float) – width of kernel density estimation (for the probability density) squared.

• Ftot_den_limit (float, optional) – Probability density limit below which data will be set to zero (this is done for numerical stability reasons. For default Ftot_den_limit, numerical difference is negligable). Defaults to 1E-10.

Returns

[pb_t, Fpbt, Fbias_window]

pb_t (array of shape (nbins,)): Probability density of window

Fpbt (array of shape (nbins,)): Force component associated with the probability density of the time-window.

Fbias_window (array of shape (nbins,)): Force component associated with the metadynamics hill deposited at the beginning of the time-window.

Return type

list

pyMFI.MFI1D.zero_to_nan(input_array)

Turns all zero elements of an array into numpy NaN (Not a Number)

Parameters

input_array (array of size (x,)) – Input array.

Returns

Input array with zero elements turned to Nan.

Return type

array

pyMFI.run_plumed module

pyMFI.run_plumed.find_alanine_dipeptide_input(initial_position_x=0.0, initial_position_y=0.0, file_extension='')

Prepares the input file for a set of initial positions. Requires traj_comp.xtc file or similar (trajecory file of a simulation that already explored specified initial positions). 1 step: Find structures that are +- 0.5nm away from intial position. 2 step: Randomly choose one of the structures. 3 step: Produce input file for alanine dipeptide simulation.

Parameters

• initial_position_x (float, optional) – x-value of initial position of simulation. Defaults to 0.0.

• initial_position_y (float, optional) – y-value of initial position of simulation. Defaults to 0.0.

• file_extension (str, optional) – Adds an extension the the structure.gro and input.tpr file. E.g. file_extension=”_1” -> structure_file=”structure_1.gro”. Defaults to “”.. Defaults to “”.

pyMFI.run_plumed.make_external_bias_1D(grid_old, FES, Ftot, grid_min_plumed=None, grid_max_plumed=None, nbins_plumed=None, file_name_extension='', return_array=None)

pyMFI.run_plumed.run_2D_Invernizzi(simulation_steps=10, sigma=0.1, height=0.5, biasfactor=10, initial_position_x=0.0, initial_position_y=0.0, hp_centre_x=0.0, hp_centre_y=0.0, hp_kappa_x=0, hp_kappa_y=0, lw_centre_x=0.0, lw_centre_y=0.0, lw_kappa_x=0, lw_kappa_y=0, uw_centre_x=0.0, uw_centre_y=0.0, uw_kappa_x=0, uw_kappa_y=0, gaus_pace=500, position_pace=0, file_extension='', external_bias_file='')

Function to run a langevin simulation (in 2D) on the Invernizzi potential. Analytical form is approx.: 1.35*x^4+1.90*x^3*y+3.93*x^2*y^2-6.44*x^2-1.90*x*y^3+5.59*x*y+1.33*x+1.35*y^4-5.56*y^2+0.90*y+18.56.

Parameters

• simulation_steps (int, optional) – Number of steps in simulation. Defaults to 100000.

• sigma (float, optional) – Gaussian width (sigma) in x-direction and y-direction of metadynamics bias. Defaults to 0.1.

• height (float, optional) – aussian height of metadynamics bias. Defaults to 0.5.

• biasfactor (int, optional) – Biasfactor of metadynamics bias. Defaults to 10.

• initial_position_x (float, optional) – x-value of initial position of simulation. Defaults to 0.0.

• initial_position_y (float, optional) – y-value of initial position of simulation. Defaults to 0.0.

• hp_centre_x (float, optional) – x-position of harmonic potential. Defaults to 0.0.

• hp_centre_y (float, optional) – y-position of harmonic potential. Defaults to 0.0.

• hp_kappa_x (int, optional) – x-force_constant of harmonic potential. Defaults to 0.

• hp_kappa_y (int, optional) – y-force_constant of harmonic potential. Defaults to 0.

• lw_centre_x (float, optional) – x-position of lower wall potential. Defaults to 0.0.

• lw_centre_y (float, optional) – y-position of lower wall potential. Defaults to 0.0.

• lw_kappa_x (int, optional) – x-force_constant of lower wall potential. Defaults to 0.

• lw_kappa_y (int, optional) – y-force_constant of lower wall potential. Defaults to 0.

• uw_centre_x (float, optional) – x-position of upper wall potential. Defaults to 0.0.

• uw_centre_y (float, optional) – y-position of upper wall potential. Defaults to 0.0.

• uw_kappa_x (int, optional) – x-force_constant of upper wall potential. Defaults to 0.

• uw_kappa_y (int, optional) – y-force_constant of upper wall potential. Defaults to 0.

• gaus_pace (int, optional) – Pace of deposition of metadynamics hills. Defaults to 500.

• position_pace (int, optional) – Pace of recording the CV in the position file. When position_pace=0, position_pace = gaus_pace/10. Defaults to 0.

• file_extension (str, optional) – Adds an extension the the position and HILLS file. E.g. file_extension=”_1” -> position_file=”position_1”. Defaults to “”.

pyMFI.run_plumed.run_alanine_dipeptide(simulation_steps, temperature=2.49, grid_min_x='-pi', grid_max_x='pi', grid_min_y='-pi', grid_max_y='pi', grid_bin_x=201, grid_bin_y=201, gaus_width_x=0.1, gaus_width_y=0.1, gaus_height=1, biasfactor=10, gaus_pace=100, position_pace=0, hp_centre_x=0.0, hp_centre_y=0.0, hp_kappa_x=0, hp_kappa_y=0, lw_centre_x=0.0, lw_centre_y=0.0, lw_kappa_x=0, lw_kappa_y=0, uw_centre_x=0.0, uw_centre_y=0.0, uw_kappa_x=0, uw_kappa_y=0, print_bias=0, file_extension='')

Function to run molecular simulation on alanine dipeptide. Requires a reference.pdb and input.tpr file.

Parameters

• simulation_steps (int) – Number of steps in simulation

• temperature (float, optional) – Temperature of simulation (units in kT). Defaults to 2.49.

• grid_min_x (str, optional) – phi-value of minimum value of grid where the bias is stored. Defaults to “-pi”.

• grid_max_x (str, optional) – phi-value of maximum value of grid where the bias is stored. Defaults to “pi”.

• grid_min_y (str, optional) – psi-value of minimum value of grid where the bias is stored. Defaults to “-pi”.

• grid_max_y (str, optional) – psi-value of maximum value of grid where the bias is stored. Defaults to “pi”.

• grid_bin_x (int, optional) – Number of distinct bins in grid on phi-axis. Defaults to 201.

• grid_bin_y (int, optional) – Number of distinct bins in grid on psi-axis. Defaults to 201.

• gaus_width_x (float, optional) – Gaussian width in phi-direction (sigma) of metadynamics bias. Defaults to 0.1.

• gaus_width_y (float, optional) – Gaussian width in psi-direction (sigma) of metadynamics bias. Defaults to 0.1.

• gaus_height (int, optional) – Gaussian height of metadynamics bias. Defaults to 1.

• biasfactor (int, optional) – Biasfactor of metadynamics bias. Defaults to 10.

• gaus_pace (int, optional) – Pace of deposition of metadynamics hills. Defaults to 100.

• position_pace (int, optional) – Pace of recording the CV in the position file. When position_pace=0, position_pace = gaus_pace/10. Defaults to 0.

• hp_centre_x (float, optional) – phi-position of harmonic potential. Defaults to 0.0.

• hp_centre_y (float, optional) – psi-position of harmonic potential. Defaults to 0.0.

• hp_kappa_x (int, optional) – phi-force_constant of harmonic potential. Defaults to 0.

• hp_kappa_y (int, optional) – psi-force_constant of harmonic potential. Defaults to 0.

• lw_centre_x (float, optional) – phi-position of lower wall potential. Defaults to 0.0.

• lw_centre_y (float, optional) – psi-position of lower wall potential. Defaults to 0.0.

• lw_kappa_x (int, optional) – phi-force_constant of lower wall potential. Defaults to 0.

• lw_kappa_y (int, optional) – psi-force_constant of lower wall potential. Defaults to 0.

• uw_centre_x (float, optional) – phi-position of upper wall potential. Defaults to 0.0.

• uw_centre_y (float, optional) – psi-position of upper wall potential. Defaults to 0.0.

• uw_kappa_x (int, optional) – phi-force_constant of upper wall potential. Defaults to 0.

• uw_kappa_y (int, optional) – psi-force_constant of upper wall potential. Defaults to 0.

• print_bias (int, optional) – When print_bias=1, the experienced bias potential and the bias force squared is printed every 100 steps in a file called “restraint”. Defaults to 0.

• file_extension (str, optional) – Adds an extension the the position and HILLS file. E.g. file_extension=”_1” -> position_file=”position_1”. Defaults to “”.

pyMFI.run_plumed.run_alanine_dipeptide_1D(simulation_steps, tors_angle='phi', temperature=2.49, grid_min='-pi', grid_max='pi', grid_bin=301, gaus_width=0.1, gaus_height=1, biasfactor=10, gaus_pace=100, position_pace=0, hp_centre=0.0, hp_kappa=0, lw_centre=0.0, lw_kappa=0, uw_centre=0.0, uw_kappa=0, print_bias=0, file_extension='')

Function to run molecular simulation on alanine dipeptide. Requires a reference.pdb and input.tpr file.

Parameters

• simulation_steps (int) – Number of steps in simulation

• temperature (float, optional) – Temperature of simulation (units in kT). Defaults to 2.49.

• grid_min_x (str, optional) – phi-value of minimum value of grid where the bias is stored. Defaults to “-pi”.

• grid_max_x (str, optional) – phi-value of maximum value of grid where the bias is stored. Defaults to “pi”.

• grid_min_y (str, optional) – psi-value of minimum value of grid where the bias is stored. Defaults to “-pi”.

• grid_max_y (str, optional) – psi-value of maximum value of grid where the bias is stored. Defaults to “pi”.

• grid_bin_x (int, optional) – Number of distinct bins in grid on phi-axis. Defaults to 201.

• grid_bin_y (int, optional) – Number of distinct bins in grid on psi-axis. Defaults to 201.

• gaus_width_x (float, optional) – Gaussian width in phi-direction (sigma) of metadynamics bias. Defaults to 0.1.

• gaus_width_y (float, optional) – Gaussian width in psi-direction (sigma) of metadynamics bias. Defaults to 0.1.

• gaus_height (int, optional) – Gaussian height of metadynamics bias. Defaults to 1.

• biasfactor (int, optional) – Biasfactor of metadynamics bias. Defaults to 10.

• gaus_pace (int, optional) – Pace of deposition of metadynamics hills. Defaults to 100.

• position_pace (int, optional) – Pace of recording the CV in the position file. When position_pace=0, position_pace = gaus_pace/10. Defaults to 0.

• hp_centre_x (float, optional) – phi-position of harmonic potential. Defaults to 0.0.

• hp_centre_y (float, optional) – psi-position of harmonic potential. Defaults to 0.0.

• hp_kappa_x (int, optional) – phi-force_constant of harmonic potential. Defaults to 0.

• hp_kappa_y (int, optional) – psi-force_constant of harmonic potential. Defaults to 0.

• lw_centre_x (float, optional) – phi-position of lower wall potential. Defaults to 0.0.

• lw_centre_y (float, optional) – psi-position of lower wall potential. Defaults to 0.0.

• lw_kappa_x (int, optional) – phi-force_constant of lower wall potential. Defaults to 0.

• lw_kappa_y (int, optional) – psi-force_constant of lower wall potential. Defaults to 0.

• uw_centre_x (float, optional) – phi-position of upper wall potential. Defaults to 0.0.

• uw_centre_y (float, optional) – psi-position of upper wall potential. Defaults to 0.0.

• uw_kappa_x (int, optional) – phi-force_constant of upper wall potential. Defaults to 0.

• uw_kappa_y (int, optional) – psi-force_constant of upper wall potential. Defaults to 0.

• print_bias (int, optional) – When print_bias=1, the experienced bias potential and the bias force squared is printed every 100 steps in a file called “restraint”. Defaults to 0.

• file_extension (str, optional) – Adds an extension the the position and HILLS file. E.g. file_extension=”_1” -> position_file=”position_1”. Defaults to “”.

pyMFI.run_plumed.run_langevin1D(simulation_steps, analytical_function='7*x^4-23*x^2', periodic='NO', initial_position=0.0, temperature=1, time_step=0.005, grid_min=- 3.0, grid_max=3.0, grid_bin=200, gaus_width=0.1, gaus_height=1, biasfactor=10, gaus_pace=100, position_pace=0, hp_centre=0.0, hp_kappa=0, lw_centre=0.0, lw_kappa=0, uw_centre=0.0, uw_kappa=0, external_bias_file='', start_sim=True)

Function to run a langevin simulation in 1 dimension. Default analytical potential: y = 7*x^4-23*x^2.

Parameters

• simulation_steps (int) – Number of steps in simulation

• analytical_function (str, optional) – The analytical function to be analysed. Defaults to “7*x^4-23*x^2”.

• periodic (str, optional) – Information wheather boundary conditions are periodic (“ON”) or not (“NO”). Defaults to “NO”.

• initial_position (float, optional) – Initial position of simulation. Defaults to 0.0.

• temperature (int, optional) – Temperature of simulation (units in kT). Defaults to 1.

• time_step (float, optional) – Length of one time step (units in ps). Defaults to 0.005.

• grid_min (float, optional) – Minimum value of grid where the bias is stored. Defaults to -3.0.

• grid_max (float, optional) – Maximum value of grid where the bias is stored. Defaults to 3.0.

• grid_bin (int, optional) – Number of distinct bins in grid. Defaults to 200.

• gaus_width (float, optional) – Gaussian width (sigma) of metadynamics bias. Defaults to 0.1.

• gaus_height (int, optional) – Gaussian height of metadynamics bias. Defaults to 1.

• biasfactor (int, optional) – Biasfactor of metadynamics bias. Defaults to 10.

• gaus_pace (int, optional) – Pace of deposition of metadynamics hills. Defaults to 100.

• position_pace (int, optional) – Pace of recording the CV in the position file. When position_pace=0, position_pace = gaus_pace/10. Defaults to 0.

• hp_centre (float, optional) – position of harmonic potential. Defaults to 0.0.

• hp_kappa (int, optional) – force_constant of harmonic potential. Defaults to 0.

• lw_centre (float, optional) – position of lower wall potential. Defaults to 0.0.

• lw_kappa (int, optional) – force_constant of lower wall potential. Defaults to 0.

• uw_centre (float, optional) – position of upper wall potential. Defaults to 0.0.

• uw_kappa (int, optional) – force_constant of upper wall potential. Defaults to 0.

• external_bias_file (str, optional) – File name or file path of external bias. Default is “”, so no file will be loaded.

pyMFI.run_plumed.run_langevin1D_plumed_fes(simulation_steps, analytical_function='7*x^4-23*x^2', periodic='NO', temperature=1, time_step=0.005, initial_position=- 1.0, gaus_width=0.1, gaus_height=1, biasfactor=10, fes_stride=0, grid_min_plumed=- 3, grid_max_plumed=3, grid_bin_plumed=301, grid_min_out=None, grid_max_out=None, grid_bin_out=None)

Function to run a langevin simulation in 1 dimension on analytical potential: y = 7*x^4-23*x^2, while also calculating the FES through plumed.

Parameters

• length (int) – Number of steps in simulation

• sigma (float, optional) – Gaussian width of metadynamics bias. Defaults to 0.1.

• height (float, optional) – Gaussian height of metadynamics bias. Defaults to 0.1.

• biasfactor (int, optional) – Biasfactor of metadynamics bias. Defaults to 10.

• fes_stride (int, optional) – Number of times the intermediate fes is calculated. Defaults to 0.

• grid_min (int, optional) – Minimum value of grid where the bias is stored. Defaults to -3.

• grid_max (int, optional) – Maximum value of grid where the bias is stored. Defaults to 3.

• grid_bin (int, optional) – Number of distinct bins in grid. Defaults to 301.

pyMFI.run_plumed.run_langevin2D(simulation_steps, analytical_function='7*x^4-23*x^2+7*y^4-23*y^2', periodic_f='NO', initial_position_x=0.0, initial_position_y=0.0, temperature=1, time_step=0.005, grid_min_x=- 3.0, grid_max_x=3.0, grid_min_y=- 3.0, grid_max_y=3.0, grid_bin_x=200, grid_bin_y=200, gaus_width_x=0.1, gaus_width_y=0.1, gaus_height=1, biasfactor=10, gaus_pace=100, hp_centre_x=0.0, hp_centre_y=0.0, hp_kappa_x=0, hp_kappa_y=0, lw_centre_x=0.0, lw_centre_y=0.0, lw_kappa_x=0, lw_kappa_y=0, uw_centre_x=0.0, uw_centre_y=0.0, uw_kappa_x=0, uw_kappa_y=0, position_pace=0, file_extension='')

Function to run a langevin simulation in 2 dimension. Default analytical potential: z = 7*x^4-23*x^2+7*y^4-23*y^2.

Parameters

• simulation_steps (int) – Number of steps in simulation

• analytical_function (str, optional) – The analytical function to be analysed. Defaults to “7*x^4-23*x^2+7*y^4-23*y^2”.

• periodic_f (str, optional) – Information wheather boundary conditions are periodic (“ON”) or not (“NO”). Defaults to “NO”.

• initial_position_x (float, optional) – x-value of initial position of simulation. Defaults to 0.0.

• initial_position_y (float, optional) – y-value of initial position of simulation. Defaults to 0.0.

• temperature (int, optional) – Temperature of simulation (units in kT). Defaults to 1.

• time_step (float, optional) – Length of one time step (units in ps). Defaults to 0.005.

• grid_min_x (float, optional) – x-value of minimum value of grid where the bias is stored. Defaults to -3.0.

• grid_max_x (float, optional) – x-value of maximum value of grid where the bias is stored. Defaults to 3.0.

• grid_min_y (float, optional) – y-value of minimum value of grid where the bias is stored. Defaults to -3.0.

• grid_max_y (float, optional) – y-value of maximum value of grid where the bias is stored. Defaults to 3.0.

• grid_bin_x (int, optional) – Number of distinct bins in grid on x-axis. Defaults to 200.

• grid_bin_y (int, optional) – Number of distinct bins in grid on y-axis. Defaults to 200.

• gaus_width_x (float, optional) – Gaussian width in x-direction (sigma) of metadynamics bias. Defaults to 0.1.

• gaus_width_y (float, optional) – Gaussian width in y-direction (sigma) of metadynamics bias. Defaults to 0.1.

• gaus_height (int, optional) – Gaussian height of metadynamics bias. Defaults to 1.

• biasfactor (int, optional) – Biasfactor of metadynamics bias. Defaults to 10.

• gaus_pace (int, optional) – Pace of deposition of metadynamics hills. Defaults to 100.

• hp_centre_x (float, optional) – x-position of harmonic potential. Defaults to 0.0.

• hp_centre_y (float, optional) – y-position of harmonic potential. Defaults to 0.0.

• hp_kappa_x (int, optional) – x-force_constant of harmonic potential. Defaults to 0.

• hp_kappa_y (int, optional) – y-force_constant of harmonic potential. Defaults to 0.

• lw_centre_x (float, optional) – x-position of lower wall potential. Defaults to 0.0.

• lw_centre_y (float, optional) – y-position of lower wall potential. Defaults to 0.0.

• lw_kappa_x (int, optional) – x-force_constant of lower wall potential. Defaults to 0.

• lw_kappa_y (int, optional) – y-force_constant of lower wall potential. Defaults to 0.

• uw_centre_x (float, optional) – x-position of upper wall potential. Defaults to 0.0.

• uw_centre_y (float, optional) – y-position of upper wall potential. Defaults to 0.0.

• uw_kappa_x (int, optional) – x-force_constant of upper wall potential. Defaults to 0.

• uw_kappa_y (int, optional) – y-force_constant of upper wall potential. Defaults to 0.

• position_pace (int, optional) – Pace of recording the CV in the position file. When position_pace=0, position_pace = gaus_pace/10. Defaults to 0.

• file_extension (str, optional) – Adds an extension the the position and HILLS file. E.g. file_extension=”_1” -> position_file=”position_1”. Defaults to “”.

pyMFI.MFI_functions module

Module contents