r"""Calculate Bending-Stretching-Twisting matrix

        It calculate elastic matrix and modulus matrix.

        .. math::

            \text{modulus matrix} = 4.1419464 \times \begin{bmatrix}
            K_{Bx}       & K_{Bx,By} & K_{Bx,S} & K_{Bx,T} \\
            K_{Bx,By}    & K_{By}    & K_{By,S} & K_{By,T} \\
            K_{Bx,S}     & K_{By,S}  & K_{S}    & K_{S,T} \\
            K_{Bx,T}     & K_{Bx,T}  & K_{S,T}  & K_{T}
            \end{bmatrix} \times L_0


        .. currentmodule:: dnaMD

        Parameters
        ----------
        bp : list
            List of two base-steps forming the DNA segment.
            For example: with ``bp=[5, 50]``, 5-50 base-step segment will be considered.

        frames : list
            List of two trajectory frames between which parameters will be extracted. It can be used to select portions
            of the trajectory. For example, with ``frames=[100, 1000]``, 100th to 1000th frame of the trajectory will be
            considered.

        paxis : str
            Axis parallel to global helical-axis(``'X'``, or ``'Y'`` or ``'Z'``). Only require when bending motions are
            included in the calculation.

        masked : bool
            ``Default=False``. To skip specific frames/snapshots.
            ``DNA.mask`` array should be set to use this functionality.
            This array contains boolean (either ``True`` or ``False``) value
            for each frame to mask the frames. Presently, mask array is
            automatically generated during :meth:`dnaMD.DNA.generate_smooth_axis` to
            skip those frames where 3D fitting curve was not successful within
            the given criteria.

        matrix : bool
            If it is ``True``, elastic constant matrix will be returned. Otherwise, by default modulus matrix will be
            returned.

        Return
        ------
        mean : numpy.ndarray
            Value of bending angles, contour length and twist angle (as 1D array) at which energy is zero. Minimum point
            on free energy landscape.

            .. math::
                \begin{bmatrix}
                    \theta^{x}_0       & \theta^{y}_0 & L_0 & \phi_0
                \end{bmatrix}

        result : numpy.ndarray
            Either elastic matrix or modulus matrix depending on ``matrix`` value.

def getStretchTwistBendModulus(self, bp, frames=None, paxis='Z', masked=True, matrix=False):
        r"""Calculate Bending-Stretching-Twisting matrix

        It calculate elastic matrix and modulus matrix.

        .. math::

            \text{modulus matrix} = 4.1419464 \times \begin{bmatrix}
            K_{Bx}       & K_{Bx,By} & K_{Bx,S} & K_{Bx,T} \\
            K_{Bx,By}    & K_{By}    & K_{By,S} & K_{By,T} \\
            K_{Bx,S}     & K_{By,S}  & K_{S}    & K_{S,T} \\
            K_{Bx,T}     & K_{Bx,T}  & K_{S,T}  & K_{T}
            \end{bmatrix} \times L_0


        .. currentmodule:: dnaMD

        Parameters
        ----------
        bp : list
            List of two base-steps forming the DNA segment.
            For example: with ``bp=[5, 50]``, 5-50 base-step segment will be considered.

        frames : list
            List of two trajectory frames between which parameters will be extracted. It can be used to select portions
            of the trajectory. For example, with ``frames=[100, 1000]``, 100th to 1000th frame of the trajectory will be
            considered.

        paxis : str
            Axis parallel to global helical-axis(``'X'``, or ``'Y'`` or ``'Z'``). Only require when bending motions are
            included in the calculation.

        masked : bool
            ``Default=False``. To skip specific frames/snapshots.
            ``DNA.mask`` array should be set to use this functionality.
            This array contains boolean (either ``True`` or ``False``) value
            for each frame to mask the frames. Presently, mask array is
            automatically generated during :meth:`dnaMD.DNA.generate_smooth_axis` to
            skip those frames where 3D fitting curve was not successful within
            the given criteria.

        matrix : bool
            If it is ``True``, elastic constant matrix will be returned. Otherwise, by default modulus matrix will be
            returned.

        Return
        ------
        mean : numpy.ndarray
            Value of bending angles, contour length and twist angle (as 1D array) at which energy is zero. Minimum point
            on free energy landscape.

            .. math::
                \begin{bmatrix}
                    \theta^{x}_0       & \theta^{y}_0 & L_0 & \phi_0
                \end{bmatrix}

        result : numpy.ndarray
            Either elastic matrix or modulus matrix depending on ``matrix`` value.
        """


        if self.esType == 'ST':
            raise KeyError(' Use dnaEY.getStretchTwistModulus for Stretching-Twisting modulus.')

        frames = self._validateFrames(frames)

        name = '{0}-{1}-{2}-{3}'.format(bp[0], bp[1], frames[0], frames[1])

        if name not in self.esMatrix:
            time, array = self.extractGlobalParameters(self.dna, bp, frames=frames, paxis=paxis, masked=masked)
            mean = np.mean(array, axis=1)
            esMatrix = np.asarray(self.getElasticMatrix(array))
            self.esMatrix[name] = esMatrix
            self.minimumPoint[name] = mean
        else:
            esMatrix = self.esMatrix[name]
            mean = self.minimumPoint[name]

        if not matrix:
            result = 4.1419464 * np.array(esMatrix) * mean[2]       # Calculate modulus
        else:
            result = esMatrix

        return mean, result