To calculate deviation in the given parameters of a Subject DNA to Reference DNA along the given axis.

    .. note:: Deviation = Reference_DNA(parameter) - Subject_DNA(parameter)

    .. warning::
            To calculate errors by using ``error = 'acf'`` or ``error = 'block'``,
            GROMACS tool ``g_analyze`` or ``gmx analyze`` should be present in ``$PATH``.

    Parameters
    ----------
    dnaRef  : :class:`DNA`
        Reference DNA
    dnaSubj : :class:`DNA`
        Subject DNA. Number of base-pairs in Reference and Subject DNA **should be** same.
    parameter : str
        Name of a base-pair or base-step or helical base-step parameter
        For details about accepted keywords, see ``parameter`` in the method :meth:`DNA.get_parameters`.

    bp : 1D list or array
        base-pairs to analyze
        Example: ::

            bp = [6]                                # bp_range = False
            bp = [4,15]                             # bp_range = True
            bp = range(4,15)                        # bp_range = False
            bp = np.arange(4,15)                    # bp_range = False
            bp = [2,5,6,7,9,12,18]                  # bp_range = False

    bp_range : bool
        ``Default=True``: As shown above, if ``True``, bp is taken as a range otherwise list or numpy array.

    axis : str
        Axis along which DNA axis is parallel. Keywords: ``X``, ``Y`` and ``Z``.

    windows  : int
        Number of bins along the axis

    err_type : str
        Method of error estimation.
        Currently accepted method as follows:

            * ``error = 'std'``   : Standard Deviation
            * ``error = 'acf'``   : Standard error using autocorrelation time (requires: ``g_analyze`` or ``gmx analyze``)
            * ``error = 'block'`` : Standard error using block averaging method (requires: ``g_analyze`` or ``gmx analyze``)

    tool : str
        Gromacs tool ``g_analyze`` or ``gmx analyze`` or ``gmx_mpi analyze`` etc.
        It will be used to calculate autocorrelation time or  block averaging error.
        It should be present in ``$PATH``

    Returns
    -------
    deviation : 1D array
        length of no. of windows; Deviation in the parameter for two given DNA
    deviation_error : 1D array
        length of no. of windows; Standard error in deviation fo each window/bin
    axis  : 1D array
        length of no. of windows; average position of window/bin along given axis
    axis_error  : 1D array
        length of no. of windows; Standard error in average position of window/bin along given axis

def dev_parameters_vs_axis(dnaRef, dnaSubj, parameter, bp, axis='Z', bp_range=True, windows=10, err_type='block', tool='gmx analyze'):
    """To calculate deviation in the given parameters of a Subject DNA to Reference DNA along the given axis.

    .. note:: Deviation = Reference_DNA(parameter) - Subject_DNA(parameter)

    .. warning::
            To calculate errors by using ``error = 'acf'`` or ``error = 'block'``,
            GROMACS tool ``g_analyze`` or ``gmx analyze`` should be present in ``$PATH``.

    Parameters
    ----------
    dnaRef  : :class:`DNA`
        Reference DNA
    dnaSubj : :class:`DNA`
        Subject DNA. Number of base-pairs in Reference and Subject DNA **should be** same.
    parameter : str
        Name of a base-pair or base-step or helical base-step parameter
        For details about accepted keywords, see ``parameter`` in the method :meth:`DNA.get_parameters`.

    bp : 1D list or array
        base-pairs to analyze
        Example: ::

            bp = [6]                                # bp_range = False
            bp = [4,15]                             # bp_range = True
            bp = range(4,15)                        # bp_range = False
            bp = np.arange(4,15)                    # bp_range = False
            bp = [2,5,6,7,9,12,18]                  # bp_range = False

    bp_range : bool
        ``Default=True``: As shown above, if ``True``, bp is taken as a range otherwise list or numpy array.

    axis : str
        Axis along which DNA axis is parallel. Keywords: ``X``, ``Y`` and ``Z``.

    windows  : int
        Number of bins along the axis

    err_type : str
        Method of error estimation.
        Currently accepted method as follows:

            * ``error = 'std'``   : Standard Deviation
            * ``error = 'acf'``   : Standard error using autocorrelation time (requires: ``g_analyze`` or ``gmx analyze``)
            * ``error = 'block'`` : Standard error using block averaging method (requires: ``g_analyze`` or ``gmx analyze``)

    tool : str
        Gromacs tool ``g_analyze`` or ``gmx analyze`` or ``gmx_mpi analyze`` etc.
        It will be used to calculate autocorrelation time or  block averaging error.
        It should be present in ``$PATH``

    Returns
    -------
    deviation : 1D array
        length of no. of windows; Deviation in the parameter for two given DNA
    deviation_error : 1D array
        length of no. of windows; Standard error in deviation fo each window/bin
    axis  : 1D array
        length of no. of windows; average position of window/bin along given axis
    axis_error  : 1D array
        length of no. of windows; Standard error in average position of window/bin along given axis

    """
    RefParam, ref_bp_idx = dnaRef.get_parameters(parameter, bp, bp_range)
    RefAxis, dummy = dnaRef.get_parameters(
        'Helical {0}-axis' .format(axis), bp, bp_range)
    SubjParam, subj_bp_idx = dnaSubj.get_parameters(parameter, bp, bp_range)

    mean_axis = np.mean(RefAxis, axis=1)
    meanRefParam = np.mean(RefParam, axis=1)
    meanSubjParam = np.mean(SubjParam, axis=1)

    maxAxis = np.amax(mean_axis)
    minAxis = np.amin(mean_axis)
    axis_range = (maxAxis - minAxis) / windows

    Ref_param_error = get_error(dnaRef.time, RefParam, len(ref_bp_idx), err_type=err_type, tool=tool)
    Ref_axis_error = get_error(dnaRef.time, RefAxis, len(ref_bp_idx), err_type=err_type, tool=tool)
    subj_param_error = get_error(dnaSubj.time, SubjParam, len(subj_bp_idx), err_type=err_type, tool=tool)

    merged_ref_param = []
    merged_subj_Param = []
    merged_Ref_param_error = []
    merged_Ref_axis_error = []
    merged_subj_param_error = []

    final_axis = []

    for i in range(windows):
        start = minAxis + (i * axis_range)
        end = start + axis_range
        idx = []
        for j in range(len(mean_axis)):
            if((start <= mean_axis[j]) and (end > mean_axis[j])):
                idx.append(j)
        if(len(idx) > 0):
            merged_ref_param.append(meanRefParam[idx])
            merged_subj_Param.append(meanSubjParam[idx])
            final_axis.append(start + (end - start) / 2)
            merged_Ref_param_error.append(Ref_param_error[idx])
            merged_Ref_axis_error.append(Ref_axis_error[idx])
            merged_subj_param_error.append(subj_param_error[idx])

    final_ref_param = []
    final_subj_param = []
    final_ref_param_error = []
    final_ref_axis_error = []
    final_subj_param_error = []
    for i in range(len(merged_ref_param)):
        final_ref_param.append(np.sum(merged_ref_param[i]))
        final_subj_param.append(np.sum(merged_subj_Param[i]))

        final_ref_axis_error.append(
            np.sqrt((merged_Ref_axis_error[i]**2).sum()))
        final_ref_param_error.append(
            np.sqrt((merged_Ref_param_error[i]**2).sum()))
        final_subj_param_error.append(
            np.sqrt((merged_subj_param_error[i]**2).sum()))

    deviation, error = get_deviation(
        final_ref_param, final_ref_param_error, final_subj_param, final_subj_param_error)

    return deviation, error, final_axis, final_ref_axis_error