r"""Deformation energy of the input DNA using local elastic properties

        The deformation energy of a base-step/s for probe DNA object with reference to
        the same base-step/s DNA present in the current DNA object.

        The deformation free energy is calculated using elastic matrix as follows

        .. math::

            G = \frac{1}{2}\mathbf{xKx^T}

        When ``helical='False'``

        .. math::
            \mathbf{K} = \mathbf{K}_{base-step}

        .. math::

            \mathbf{x} =  \begin{bmatrix}
                              (Dx_{i}-Dx_0)  &  (Dy_i - Dy_0) & (Dz_i - Dz_0) & (\tau_i - \tau_0) &
                              (\rho_i - \rho_0) & (\omega_i - \omega_0)
                          \end{bmatrix}


        When ``helical='True'``

        .. math::
            \mathbf{K} = \mathbf{K}_{helical-base-step}

        .. math::
            \mathbf{x} =  \begin{bmatrix}
                              (dx_{i}-dx_0)  &  (dy_i - dy_0) & (h_i - h_0) & (\eta_i - \eta_0) &
                              (\theta_i - \theta_0) & (\Omega_i - \Omega_0)
                          \end{bmatrix}

        .. 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.

        complexDna : :class:`dnaMD.DNA`
            Input :class:`dnaMD.DNA` instance for which deformation energy will be calculated.

        freeDnaFrames : list
            To select a trajectory segment of current (free) DNA data.
            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.

        boundDnaFrames : list
            To select a trajectory segment of input (bound) DNA data.
            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.

        helical : bool
            If ``helical=True``, elastic matrix for **helical base-step** parameters are calculated. Otherwise,
            by default, elastic matrix for **base-step** parameters are calculated.

        unit : str
            Unit of energy. Allowed units are: ``'kT', 'kJ/mol' and 'kcal/mol'``.

        which : str or list
            For which motions (degrees of freedom), energy should be calculated. It should be either a list containing
            terms listed below or"all" for all energy terms.

            Following keywords are available:
                * ``'full'`` : Use entire elastic matrix -- all parameters with their coupling
                * ``'diag'`` : Use diagonal of elastic matrix -- all motions but no coupling
                * ``'shift'`` or ``'x-disp'``
                * ``'slide'`` or ``'y-idsp'``
                * ``'rise'`` or ``'h-rise'``
                * ``'tilt'`` or ``'inclination'``
                * ``'roll'`` or ``'tip'``
                * ``'twist'`` or ``'h-twist'``

        outFile : str
            Output file in csv format.

        Returns
        -------
        time : numpy.ndarray
            1D array containing time values.

        energy : dict of numpy.ndarray
            Dictionary of 1D array of shape (nframes) containing energy terms requested for DNA.

def getLocalDeformationEnergy(self, bp, complexDna, freeDnaFrames=None, boundDnaFrames=None, helical=False,
                                  unit='kT', which='all', outFile=None):
        r"""Deformation energy of the input DNA using local elastic properties

        The deformation energy of a base-step/s for probe DNA object with reference to
        the same base-step/s DNA present in the current DNA object.

        The deformation free energy is calculated using elastic matrix as follows

        .. math::

            G = \frac{1}{2}\mathbf{xKx^T}

        When ``helical='False'``

        .. math::
            \mathbf{K} = \mathbf{K}_{base-step}

        .. math::

            \mathbf{x} =  \begin{bmatrix}
                              (Dx_{i}-Dx_0)  &  (Dy_i - Dy_0) & (Dz_i - Dz_0) & (\tau_i - \tau_0) &
                              (\rho_i - \rho_0) & (\omega_i - \omega_0)
                          \end{bmatrix}


        When ``helical='True'``

        .. math::
            \mathbf{K} = \mathbf{K}_{helical-base-step}

        .. math::
            \mathbf{x} =  \begin{bmatrix}
                              (dx_{i}-dx_0)  &  (dy_i - dy_0) & (h_i - h_0) & (\eta_i - \eta_0) &
                              (\theta_i - \theta_0) & (\Omega_i - \Omega_0)
                          \end{bmatrix}

        .. 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.

        complexDna : :class:`dnaMD.DNA`
            Input :class:`dnaMD.DNA` instance for which deformation energy will be calculated.

        freeDnaFrames : list
            To select a trajectory segment of current (free) DNA data.
            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.

        boundDnaFrames : list
            To select a trajectory segment of input (bound) DNA data.
            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.

        helical : bool
            If ``helical=True``, elastic matrix for **helical base-step** parameters are calculated. Otherwise,
            by default, elastic matrix for **base-step** parameters are calculated.

        unit : str
            Unit of energy. Allowed units are: ``'kT', 'kJ/mol' and 'kcal/mol'``.

        which : str or list
            For which motions (degrees of freedom), energy should be calculated. It should be either a list containing
            terms listed below or"all" for all energy terms.

            Following keywords are available:
                * ``'full'`` : Use entire elastic matrix -- all parameters with their coupling
                * ``'diag'`` : Use diagonal of elastic matrix -- all motions but no coupling
                * ``'shift'`` or ``'x-disp'``
                * ``'slide'`` or ``'y-idsp'``
                * ``'rise'`` or ``'h-rise'``
                * ``'tilt'`` or ``'inclination'``
                * ``'roll'`` or ``'tip'``
                * ``'twist'`` or ``'h-twist'``

        outFile : str
            Output file in csv format.

        Returns
        -------
        time : numpy.ndarray
            1D array containing time values.

        energy : dict of numpy.ndarray
            Dictionary of 1D array of shape (nframes) containing energy terms requested for DNA.

        """
        if helical:
            energyTerms =  ['full', 'diag', 'x-disp', 'y-disp', 'h-rise', 'inclination', 'tip', 'h-twist']
        else:
            energyTerms = ['full', 'diag', 'shift', 'slide', 'rise', 'tilt', 'roll', 'twist']

        if isinstance(which, str):
            if which != 'all':
                raise ValueError('Either use "all" or use list of terms from this {0} list \n.'.format(energyTerms))
            else:
                which = energyTerms
        elif isinstance(which, list):
            for key in which:
                if key not in energyTerms:
                    raise ValueError('{0} is not a supported keyword.\n Use from the following list: \n{1}'.format(
                        which, energyTerms))
        else:
            raise ValueError('Either use "all" or use list of terms from this {0} list \n.'.format(energyTerms))

        means, esMatrix = self.calculateLocalElasticity(bp, frames=freeDnaFrames, helical=helical, unit=unit)
        time, array = self.extractLocalParameters(complexDna, bp, frames=boundDnaFrames, helical=helical)

        # Initialize energy dictionary
        energyOut = OrderedDict()
        for key in which:
            energyOut[key] = []

        for i in range(array[0].shape[0]):
            vec = array[:, i]
            diff = vec - means
            for key in which:
                t_energy = self._calcLocalEnergy(diff, esMatrix, key)
                energyOut[key].append(t_energy)

        for key in which:
            energyOut[key] = np.asarray(energyOut[key])

        # Write output file
        if outFile is not None:
            with open(outFile, 'w') as fout:
                fout.write('#Time')
                for name in which:
                    fout.write(', {0}'.format(name))
                fout.write('\n')
                for t in range(len(time)):
                    fout.write('{0:.3f}'.format(time[t]))
                    for name in which:
                        fout.write(', {0:.5f}'.format(energyOut[name][t]))
                    fout.write('\n')

        return time, energyOut