#  Copyright (C) 2020 The authors of Py-ChemShell
#
#  This file is part of Py-ChemShell.
#
#  Py-ChemShell is free software: you can redistribute it and/or modify
#  it under the terms of the GNU Lesser General Public License as
#  published by the Free Software Foundation, either version 3 of the
#  License, or (at your option) any later version.
#
#  Py-ChemShell is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU Lesser General Public License for more details.
#
#  You should have received a copy of the GNU Lesser General Public
#  License along with Py-ChemShell.  If not, see
#  <http://www.gnu.org/licenses/>.

__author__ = 'You Lu <you.lu@stfc.ac.uk>'
__date__ = 'Feb. 2020'

import chemsh
from sys import stdout

class Polynucleotide():
    ''''''

    def __new__(cls, *args, **kwargs):
        ''''''

        from .fragbio  import FragmentBio 

        inst  = FragmentBio.__new__(FragmentBio)
        inst2 = FragmentBio.__new__(FragmentBio)

        # necessary
        inst.__init__()
        inst2.__init__()

        inst = cls._init(inst, inst2, *args, **kwargs)

        return inst


    def _init(strand1Obj,
              strand2Obj,
              strand1,
              strand2      = '',
              ff           = 'charmm',
              conformation = 'B',
              backbone     = 'drew1981',
              periodic     =   10.4,
              alpha        =  -30.0,
              beta         =  136.0,
              gamma        =   31.0,
              delta        =  143.0,
              epsilon      = -141.0,
              zeta         = -161.0,
              chi          =  -98.0):
        ''''''

        from numpy           import all, arange, arccos, arctan2, array, char, clip, cross, dot, einsum, float64
        from numpy           import intersect1d, isnan, linalg, nan, rad2deg, searchsorted, sign, stack, sum, unique, zeros
        from os              import path
        from .fragment       import Fragment
        from ..data          import charmm, forcefield, unitconvert
        from ..dl_py2f.utils import fileutils

        # Ref.: Principles of Nucleic Acid Structure
        #       Authors (view affiliations)
        #       Wolfram Saenger
        # A-DNA
        if conformation.strip().lower() == 'a':
            # Ref. 304 in Saenger's book:
            if backbone.strip().lower() in [ 'arnott1980', '' ]:
                alpha, beta, gamma, delta, epsilon, zeta, chi = -90.0, 149.0, 47.0, 83.0, -175.0, -45.0, -154.0
            # Ref. 843 in Saenger's book:
            #     S. Arnott and D. W. L. Hukiris (1972), Optimised parameters for A-DNA
            #     and B-DNA. Biochem. Biophys. Res. Commun. 47, 1504-1510. 
            elif backbone.strip().lower() == 'arnott1972':
                alpha, beta, gamma, delta, epsilon, zeta, chi = -75.0, -179.0, 59.0, 79.0, -155.0, -67.0, -158.0
        # B-DNA
        elif conformation.strip().lower() == 'b':
            # Ref. 304 in Saenger's book:
            #     S. Arnott, R. Chandrasekaran, D. L. Birdsall, A. G. W. Leslie, and R. L.
            #     Ratcliff (1980), Left handed DNA helices. Nature 283, 743-745. 
            if backbone.strip().lower() in [ 'arnott1980', '' ]:
                alpha, beta, gamma, delta, epsilon, zeta, chi = -41.0, 136.0, 38.0, 139.0, -133.0, -157.0, -102.0
            # Ref. 856 in Saenger's book:
            #     S. Amott and D. W. L. Hukins (1973), Refinement of the structure of
            #     B-DNA and implications for the analysis of X-ray diffraction data from
            #     Fibers of biopolymers. J. Mol. BioI. 81, 93-105. 
            elif backbone.strip().lower() == 'arnott1973':
                alpha, beta, gamma, delta, epsilon, zeta, chi = -46.0, -147.0, 36.0, 157.0, 155.0, -96.0, -98.0
            # H. R. Drew, R. M. Wing, T. Takano, C. Broka, S. Tanaka, K. Itakura and R. E. Dickerson,
            # Proceedings of the National Academy of Sciences, 1981, 78, 2179-2183.
            elif backbone.strip().lower() == 'drew1981':
                alpha, beta, gamma, delta, epsilon, zeta, chi = -63.0, 171.0, 54.0, 123.0, -169.0, -108.0, -117.0
        elif conformation.strip().lower() == 'c':
            # Ref. 857 in Saenger's book:
            #    S. Arnott and E. Selsing (1975), The conformation of C-DNA. J. Mol. BioI.
            #    98, 265-269. 
            if backbone.strip().lower() in [ 'arnott1975', '' ]:
                alpha, beta, gamma, delta, epsilon, zeta, chi = -39.0, -160.0, 37.0, 157.0, 161.0, -106.0, -97.0
        # TODO: Z-DNA
        elif conformation.lower() == 'z':
              pass


        def __getDeletedMatrix(_blkObj, _todelete):

            # save the Z-matrix of O3' and P and delete O3'and P
            _zmatrix_deleted = _blkObj._master['zmatrix'][_todelete]
            _blkObj.delete(_todelete, mute=True)

            # recalculate Z-matrix after deleting (with implicit reconnection)
            _blkObj._calcZMatrix = True
            _blkObj.c2z()

            return _zmatrix_deleted


        def __setPO3p(_strandObj, _block, _nshift, _range_blk, _zmatrix_deleted_hydroxy):
            ''''''

            # P-O3'
            P   = _strandObj.select(atoms='P'   , residues=_block, indices=_range_blk)
            O3p = _strandObj.select(atoms='O3\'', residues='X'   , indices=_range_blk)
            P_O3p = array([P, O3p]).ravel()

            # connect P-O3' with C3p: loop over fields (bond, angle, dihedral)
            for f in _strandObj._master['zmatrix'].dtype.fields:
                _strandObj._master['zmatrix'][f][P_O3p] = _zmatrix_deleted_hydroxy[f][:2]

            # shift the indices
            for f in ['j','k','l']:
                _strandObj._master['zmatrix'][f][P_O3p] += _nshift

            # connect P with O3'
            _strandObj._master['zmatrix']['j'][P] = O3p


        def __setDihedrals0(__strandObj, __C5p, __C4p, __C3p, __O3p, __P, __O5p):
            '''Set the backbone torsions of the first 5'-terminal block'''

            # alpha: C5'-O5'-P-O3'
            __strandObj._master['zmatrix']['dihedral'][__C5p] = alpha
            # beta: C4'-C5'-O5'-P
            __strandObj._master['zmatrix']['dihedral'][__C4p] = beta
            # gamma: C3'-C4'-C5'-O5'
            __strandObj._master['zmatrix']['dihedral'][__C3p] = gamma


        def __setDihedrals(__strandObj, __C5p, __C4p, __C3p, __O3p, __P, __O5p):
            '''Set the backbone torsions'''

            # alpha: C5'-O5'-P-O3'
            __strandObj._master['zmatrix']['dihedral'][__C5p] = alpha
            # beta: C4'-C5'-O5'-P
            __strandObj._master['zmatrix']['dihedral'][__C4p] = beta
            # gamma: C3'-C4'-C5'-O5'
            __strandObj._master['zmatrix']['dihedral'][__C3p] = gamma
            # delta: O3'-C3'-C4'-C5'
            __strandObj._master['zmatrix']['dihedral'][__O3p] = delta
            # epsilon: P-O3'-C3'-C4'
            __strandObj._master['zmatrix']['dihedral'][__P]   = epsilon
            # zeta: O5'-P-O3'-C3'
            __strandObj._master['zmatrix']['dihedral'][__O5p] = zeta
            

        def __setChi(__strandObj, __C2, __C4, __N9):
            '''Set the chi dihedral'''

            # purines (A and G): C2-N1-C1'-O4'
            if isnan(__N9):
                __strandObj._master['zmatrix']['dihedral'][__C2] = chi

            # pyrimidines (T and C): C4-N9-C1'-O4'
            else:
                __strandObj._master['zmatrix']['dihedral'][__C4] = chi
         

        def _getHydroxy(_blkObj, _block):
            '''Saved current one'''

            # 3'-terminal: C3'-O3'-P
            _hydroxy_H = _blkObj.select(atoms='P'   , residues='X')
            _hydroxy_O = _blkObj.select(atoms='O3\'', residues=_block)
            _hydroxy = array([_hydroxy_H, _hydroxy_O]).ravel()

            return _hydroxy


        def __getIndices(__strandObj, __range_blk, __block):
            '''Get atom indices on sugar'''

            O5p = __strandObj.select(atoms='O5\'', indices=__range_blk)[0]
            C5p = __strandObj.select(atoms='C5\'', indices=__range_blk, residues=__block)[0]
            C4p = __strandObj.select(atoms='C4\'', indices=__range_blk, residues='DEO1') [0]
            C3p = __strandObj.select(atoms='C3\'', indices=__range_blk, residues=__block)[0]
            O3p = __strandObj.select(atoms='O3\'', indices=__range_blk, residues='X')    [0]
            C1p = __strandObj.select(atoms='C1\'', indices=__range_blk, residues='DEO1') [0]
            C2p = __strandObj.select(atoms='C2\'', indices=__range_blk, residues='DEO1') [0]
            O4p = __strandObj.select(atoms='O4\'', indices=__range_blk, residues='DEO1') [0]
            P   = __strandObj.select(atoms='P'   , indices=__range_blk, residues=__block)[0]

            return O5p, C5p, C4p, C3p, O3p, C1p, C2p, O4p, P


        def __getMoreIndices(__strandObj, __range_blk, __block):
            '''Get atom indices on purine and pyrimidine'''

            # YL 26/02/2020: use strict=True to exlude atoms with prime
            N1  = __strandObj.select(atoms='N1', indices=__range_blk, residues=__block, strict=True)[0]
            C2  = __strandObj.select(atoms='C2', indices=__range_blk, residues=__block, strict=True)[0]
            N3  = __strandObj.select(atoms='N3', indices=__range_blk, residues=__block, strict=True)[0]
            C4  = __strandObj.select(atoms='C4', indices=__range_blk, residues=__block, strict=True)[0]
            C5  = __strandObj.select(atoms='C5', indices=__range_blk, residues=__block, strict=True)[0]
            C6  = __strandObj.select(atoms='C6', indices=__range_blk, residues=__block, strict=True)[0]
            # only A and G contain N7, C8, and N9
            try:
                N7 = __strandObj.select(atoms='N7', indices=__range_blk, residues=__block, strict=True)[0]
                C8 = __strandObj.select(atoms='C8', indices=__range_blk, residues=__block, strict=True)[0]
                N9 = __strandObj.select(atoms='N9', indices=__range_blk, residues=__block, strict=True)[0]
            except:
                N7 = C8 = N9 = nan

            return N1, C2, N3, C4, C5, C6, N7, C8, N9


        def _getBlkObj(_block):
            '''Construct the unit block of nucleotide'''

            from .fragbio import FragmentBio 
            from ..io     import file

            _blkObj = FragmentBio()
            file.parse(path.join(prefix, _block+'MP.pdb'), 'pdb', receiver=_blkObj, mute=True)
            _blkObj._calcZMatrix = True
            _blkObj.c2z()

            return _blkObj


        def _setStrandObj0(_strandObj, _blkObj, _block, _nshift):
            ''''''

            # indices of O5', C5', C4', C3', O3', C1', C2', O4', P in the current strand
            O5p, C5p, C4p, C3p, O3p, C1p, C2p, O4p, P = __getIndices(_strandObj, range(_strandObj.natoms), _block)
            N1, C2, N3, C4, C5, C6, N7, C8, N9 = __getMoreIndices(_strandObj, range(_strandObj.natoms), _block)

            # re-orient C1' to C1'-C2'-C3'-C4' and O4' to O4'-C1'-C2'-C3' to keep the sugar's 5-membered ring
            _strandObj.reorientDihedral(C1p, C2p, C3p, C4p)
            _strandObj.reorientDihedral(O4p, C1p, C2p, C3p)

            # re-orient to keep the purine (A and G) or pyrimidine (T and C) ring
            # YL 26/02/2020: we have to do this before we implement ring detecting in c2z()
            if isnan(N9):
                _strandObj.reorientDihedral(C5, C4, N3, C2)
                _strandObj.reorientDihedral(C6, C5, C4, N3)
            else:
                _strandObj.reorientDihedral(N1, C2, N3, C4)
                _strandObj.reorientDihedral(C6, N1, C2, N3)
                _strandObj.reorientDihedral(C5, C6, N1, C2)
                _strandObj.reorientDihedral(N7, C5, C6, N1)
                _strandObj.reorientDihedral(C8, N7, C5, C6)

            # dihedrals (alpha-zeta) determining the backbone
            __setDihedrals0(_strandObj, C5p, C4p, C3p, O3p, P, O5p)

            # chi
            __setChi(_strandObj, C2, C4, N9)


        def _setStrandObj(_strandObj, _blkObj, _block, _nshift, _range_blk, _zmatrix_deleted_hydroxy):
            ''''''

            # P-O3'
            __setPO3p(_strandObj, _block, _nshift, _range_blk, _zmatrix_deleted_hydroxy)

            # indices of O5', C5', C4', C3', O3', C1', C2', O4', P in the current strand
            O5p, C5p, C4p, C3p, O3p, C1p, C2p, O4p, P = __getIndices(_strandObj, _range_blk, _block)
            N1, C2, N3, C4, C5, C6, N7, C8, N9 = __getMoreIndices(_strandObj, _range_blk, _block)

            # connect to C3' of the preceding block (not the current one)
            C3p_blk = _blkObj.select(atoms='C3\'', residues=_block)
            _strandObj._master['zmatrix']['l'][O5p] = C3p_blk[-1] + _nshift

            # re-orient C1' to C1'-C2'-C3'-C4' and O4' to O4'-C1'-C2'-C3' to keep the sugar's 5-membered ring
            _strandObj.reorientDihedral(C1p, C2p, C3p, C4p)
            _strandObj.reorientDihedral(O4p, C1p, C2p, C3p)

            # re-orient to keep the purine (A and G) or pyrimidine (T and C) ring
            # YL 26/02/2020: we have to do this before we implement ring detecting in c2z()
            if isnan(N9):
                _strandObj.reorientDihedral(C5, C4, N3, C2)
                _strandObj.reorientDihedral(C6, C5, C4, N3)
            else:
                _strandObj.reorientDihedral(N1, C2, N3, C4)
                _strandObj.reorientDihedral(C6, N1, C2, N3)
                _strandObj.reorientDihedral(C5, C6, N1, C2)
                _strandObj.reorientDihedral(N7, C5, C6, N1)
                _strandObj.reorientDihedral(C8, N7, C5, C6)

            # dihedrals (alpha-zeta) determining the backbone
            __setDihedrals(_strandObj, C5p, C4p, C3p, O3p, P, O5p)

            # chi
            __setChi(_strandObj, C2, C4, N9)


        # 3'-terminal (replace with a hydrogen at an arbitrary bond length of 1 Angstrom)
        def _3p_terminal(_strandObj, inverse_z=False):
            ''''''

            P_dummy = _strandObj.select(atoms='P', residues='X')
            _strandObj.names[P_dummy] = 'H3T'
            _strandObj.znums[P_dummy] = 1
            _strandObj._master['zmatrix']['bond'][P_dummy] = 1.89

            # calculate the Cartesian coordinates
            if chemsh.CHEMSH_DEBUG > 4:
                _strandObj.printZMatrix(start=0)
            _strandObj.z2c(inverse_z=inverse_z)

            # delete all dummy C3'
            C3p_dummy = _strandObj.select(atoms='C3\'', residues='X')
            _strandObj.delete(C3p_dummy, mute=True)


        def _getAxis(_strandObj, _nblocks):
            '''Get the central axis (only valid when nblocks>=3)'''

            from numpy  import linalg, mean, multiply, subtract
            from ..base import geometry

            # let's use the O3'-O3' vectors
            vectors = zeros(shape=(_nblocks,3), dtype=float64)
            O3p_all = _strandObj.select(atoms='O3\'')
            P_all   = _strandObj.select(atoms='P')
            for iblock in range(_nblocks):
                # do not normalise
                vectors[iblock] = _strandObj.coords[O3p_all[iblock]]-_strandObj.coords[O3p_all[iblock+1]]
            
            # least-squares solution to fit an equidistant central vector
            A = (vectors[1:]-vectors[0])*2.0
            B = sum(vectors[1:]**2-vectors[0]**2, axis=1)
            axis, residuals, rank, singulars = linalg.lstsq(A, B, rcond=None)

            # normalise
            # direction: from 5' to 3'
            axis = -axis/linalg.norm(axis)

            # the strand's length by projecting the longest O3'-O3' vector onto the central axis
            vector_P_P = _strandObj.coords[P_all[0]]-_strandObj.coords[P_all[-1]]
            P_on_axis = dot(dot(axis, vector_P_P), axis)
            vector_O3p_O3p = _strandObj.coords[O3p_all[0]]-_strandObj.coords[O3p_all[-1]]
            O3p_on_axis = dot(dot(axis, vector_O3p_O3p), axis)

            # now we have got the direction using vector but we still need a point in the cartesian coordinates
            # all O3' projected on plane axis[0]*x + axis[1]*y + axis[2]*z = 0 containing the origin point [0.0,0.0,0.0]
            O3p_on_plane = geometry.projectOntoPlane(_strandObj.coords[O3p_all], a=axis[0], b=axis[1], c=axis[2])
            P_on_plane = geometry.projectOntoPlane(_strandObj.coords[P_all], a=axis[0], b=axis[1], c=axis[2])

            # least-squares solution to fit an equidistant centre point
            A = (O3p_on_plane[1:]-O3p_on_plane[0])*2.0
            B = sum(O3p_on_plane[1:]**2-O3p_on_plane[0]**2, axis=1)
            origin, residuals, rank, singulars = linalg.lstsq(A, B, rcond=None)

            # get the first and last P-origin vectors
            P_o_0 = origin - P_on_plane[0]
            P_o_1 = origin - P_on_plane[1]
            P_o_n = origin - P_on_plane[-1]

            # diameter in a.u.
            diameter = mean(linalg.norm(subtract(P_on_plane, origin), axis=1))*2.0

            return axis, origin, linalg.norm(O3p_on_axis), P_o_0, P_o_1, P_o_n, diameter


        def _getSignedAngle(v1, v2):
            '''Return the signed angle between two vectors'''

            # normalise
            v1, v2 = v1/linalg.norm(v1), v2/linalg.norm(v2)

            return rad2deg(sign(linalg.det(stack((v1[-2:],v2[-2:]))))*arccos(clip(dot(v1,v2),-1.0,1.0)))


        def _renameStrand(_strandObj, _dict_stacks, _ffObj, _segname, _chainID):
            ''''''
            from numpy import array

            # rename the segment
            _strandObj.segments = _segname

            # rename the chainID
            _strandObj.chainIDs = _chainID

            # rename the residues
            for stack, resname in _dict_stacks.items():
                try:
                    _strandObj.residues[stack] = _ffObj.nucleotide_aliases[resname]
                # the final block was deleted with the 3'-terminal and has a changed range
                except IndexError:
                    _strandObj.residues[stack[0]:] = _ffObj.nucleotide_aliases[resname]


        # path to built-in block structures
        prefix = path.join(path.dirname(__file__), '..', 'data', 'nucleotides')

        counterbases = { 'dA':'dT',
                         'dT':'dA',
                         'dG':'dC',
                         'dC':'dG',
                         'dU':'dA',
                          'A': 'U',
                          'U': 'A',
                          'G': 'C',
                          'C': 'G',
                       }

        # automatic counter-strand
        if strand2.lower() in [ 'auto', 'automatic' ]:
            strand2 = ''
        else:
            counterbases = {}

        # delete the default exisiting dummy atom
        strand1Obj.delete(0, mute=True)
        strand2Obj.delete(0, mute=True)

        strand1_lower = iter(strand1.lower())

        # for finally renaming the residues
        ffObj = forcefield.getFFObj(ff)
        stacks_strand1 = {}
        stacks_strand2 = {}
        blocks = []

        # construct strand1
        for i, n in enumerate(strand1_lower):
            if n == 'd':
                try:
                    block = n+next(strand1_lower).upper()
                except StopIteration:
                    print("\n ERROR: please check your sequence \"{}\"\n".format(strand1))
                    exit(999)
                block2 = counterbases.get(block, '')
                strand2 += block2
            else:
                block = n.upper()
                block2 = counterbases.get(block, '')
                strand2 += block2
            blocks += [block]

            # initialise block and calculate Z-matrix
            blkObj = _getBlkObj(block)

            # do the following for non-3'-terminal blocks: delete 3'-hydroxy
            if strand1_lower.__length_hint__():
                # 3'-terminal: C3'-O3'-P
                hydroxy = _getHydroxy(blkObj, block)
                # keep the deleted part for later use
                if i == 0:
                    zmatrix_deleted_hydroxy  = __getDeletedMatrix(blkObj, hydroxy)
                    nshift  = 0
                else:
                    zmatrix_deleted_hydroxy_previous, zmatrix_deleted_hydroxy  = zmatrix_deleted_hydroxy, __getDeletedMatrix(blkObj, hydroxy)
            # 3'-terminal block: use the deleted part of the previous block
            else:
                zmatrix_deleted_hydroxy_previous  = zmatrix_deleted_hydroxy

            strand1Obj = strand1Obj.merge(blkObj, reconnect=False)
            # prevent recalculating for now
            strand1Obj._calcZMatrix = False

            if i:

                range_blk = strand1Obj.stacks['FragmentBio1']
                # set resIDs (temporary for the counter-strand which will be changed later)
                strand1Obj.resIDs[range_blk] = i+1
                # YL 24/01/2022: the bridge O3' should belong to the previous block
                strand1Obj.resIDs[range_blk[0]+1] = i

                # shift the indices
                for f in ['j','k','l']:
                    # YL 09/01/2021: we only shift -1 entries except for which indices were already shifted in the above merge() indices were already shifted
                    masks_minus_one = strand1Obj._master['zmatrix'][f][range_blk] == -1
                    strand1Obj._master['zmatrix'][f][array(range_blk)[masks_minus_one]] += range_blk[0]

                _setStrandObj(strand1Obj, blkObj, block, nshift, range_blk, zmatrix_deleted_hydroxy_previous)

                # save the offset for next block
                nshift = strand1Obj.natoms-blkObj.natoms

            # skip the first block (viz. i==0)
            else:
                _setStrandObj0(strand1Obj, blkObj, block, nshift)
                strand1Obj.resIDs = 1

            # remember the positions of the blocks
            from numpy import array
            rbuff = strand1Obj.stacks['FragmentBio1']
            if i:
                stacks_strand1.update({range(array(rbuff)[0]+1-i,array(rbuff)[-1]-i+1):block})
            else:
                stacks_strand1.update({range(array(rbuff)[0]-i,array(rbuff)[-1]-i+1):block})

        # check if DNA or RNA
        if all(~char.startswith(blocks, 'd')):
            segname = 'RNA'
        elif all(char.startswith(blocks, 'd')):
            segname = 'DNA'
        else:
            segname = 'NUCL'

        # loop over strand2 in the inverse way because our code is designed to start from the 5'-end
        if strand2:
            nblocks = unique(strand1Obj.resIDs).size
        else:
            nblocks = 0

        # construct the counter-strand 
        for iblock in range(nblocks):
            # starting from 1
            block = counterbases[intersect1d(list(counterbases.keys()),
                                             char.strip(char.decode(strand1Obj.residues[strand1Obj.select(resIDs=nblocks-iblock)])))[0]]

            # initialise block and calculate Z-matrix
            blkObj = _getBlkObj(block)

            # do the following for non-3'-terminal blocks: delete 3'-hydroxy
            if iblock != nblocks-1:
                # 3'-terminal: C3'-O3'-P
                hydroxy = _getHydroxy(blkObj, block)
                # keep the deleted part for later use
                if iblock == 0:
                    zmatrix_deleted_hydroxy  = __getDeletedMatrix(blkObj, hydroxy)
                    nshift  = 0
                else:
                    zmatrix_deleted_hydroxy_previous, zmatrix_deleted_hydroxy  = zmatrix_deleted_hydroxy, __getDeletedMatrix(blkObj, hydroxy)
            # 3'-terminal block: use the deleted part of the previous block
            else:
                zmatrix_deleted_hydroxy_previous  = zmatrix_deleted_hydroxy
                
            strand2Obj = strand2Obj.merge(blkObj, reconnect=False)
            # prevent recalculating for now
            strand2Obj._calcZMatrix = False

            # skip the first block (viz. i==0)
            if iblock:

                range_blk = strand2Obj.stacks['FragmentBio1']

                # set resIDs (temporary for the counter-strand which will be changed later)
                strand2Obj.resIDs[range_blk] = nblocks+1+iblock
                # YL 24/01/2022: the bridge O3' should belong to the previous block
                strand2Obj.resIDs[range_blk[0]+1] = nblocks+iblock

                # shift the indices
                for f in ['j','k','l']:
                    # YL 09/01/2021: we only shift -1 entries except for which indices were already shifted in the above merge() indices were already shifted
                    masks_minus_one = strand2Obj._master['zmatrix'][f][range_blk] == -1
                    strand2Obj._master['zmatrix'][f][array(range_blk)[masks_minus_one]] += range_blk[0]

                _setStrandObj(strand2Obj, blkObj, block, nshift, range_blk, zmatrix_deleted_hydroxy_previous)

                # save the offset for next block
                nshift = strand2Obj.natoms-blkObj.natoms

            else:
                _setStrandObj0(strand2Obj, blkObj, block, nshift)
                strand2Obj.resIDs = nblocks+1
            # remember the positions of the blocks
            rbuff = strand2Obj.stacks['FragmentBio1']
            if iblock:
                stacks_strand2.update({range(array(rbuff)[0]+1-iblock,array(rbuff)[-1]-iblock+1):block})
            else:
                stacks_strand2.update({range(array(rbuff)[0]-iblock,array(rbuff)[-1]-iblock+1):block})

        # finish the 3'-terminal and recalculate cartesian
        _3p_terminal(strand1Obj)

        nblocks = unique(strand1Obj.resIDs).size
        # find the orientation
        if nblocks >= 3:
            axis_strand1, origin1, len_strand1, PO0_strand1, PO1_strand1, POn_strand1, diameter1 = _getAxis(strand1Obj, nblocks)

        # pair the two strands
        if strand2 and nblocks >= 3:

            _3p_terminal(strand2Obj)

            axis_strand2, origin2, len_strand2, PO0_strand2, PO1_strand2, POn_strand2, diameter2 = _getAxis(strand2Obj, nblocks)

            # align them together to share a common central axis (two opposite vectors)
            strand1Obj.translate(vector=-origin1)
            strand2Obj.translate(vector=-origin2)

            # rotate strand2 to point oppositely (origin1/origin2 are normal vectors)
            strand2Obj.rotate(theta=180.0, vector=origin1)

            # translate by a distance of the length
            strand2Obj.translate(vector=axis_strand1*(len_strand1+len_strand2)*0.5)

            # get a new origin
            axis_strand2, origin2, len_strand2, PO0_strand2, PO1_strand2, POn_strand2, diameter2 = _getAxis(strand2Obj, nblocks)

            # rotate to form H-bonds
            phase0 = _getSignedAngle(PO0_strand1, POn_strand2)
            if phase0 < 0.0:
                phase0 = 360.0 + phase0
            rotation = phase0 - 200.0
            # counter-clockwise
            strand2Obj.rotate(theta=rotation, vector=axis_strand2)

            # debug print
            if chemsh.CHEMSH_DEBUG > 2:
                axis_strand2, origin2, len_strand2, PO0_strand2, PO1_strand2, POn_strand2, diameter2 = _getAxis(strand2Obj, nblocks)
                phase1 = _getSignedAngle(PO0_strand1, POn_strand2)
                if phase1 < 0.0:
                    phase1 = 360.0 + phase1
                print("\n >>> Counter-strand rotated about the central axis by {:<.1f} degrees: {:<.1f} -> {:<.1f}".format(rotation, phase0, phase1))
                fileutils.archiveFile('_strand1.xyz')
                fileutils.archiveFile('_strand2.xyz')
                fileutils.archiveFile('_strand1.pdb')
                fileutils.archiveFile('_strand2.pdb')
                strand1Obj.save('_strand1.xyz')
                strand2Obj.save('_strand2.xyz')
                strand1Obj.save('_strand1.pdb')
                strand2Obj.save('_strand2.pdb')

            hbonds = { 'dA':['N1','H61'],
                       'dT':['H3','O4' ],
                       'dG':['O6' ,'H1','H21'],
                       'dC':['H41','N3','O2' ],
                       'A' :['U'],
                       'U' :['A'],
                       'G' :['C'],
                       'C' :['G'],
                     }

            # rename the residues
            _renameStrand(strand1Obj, stacks_strand1, ffObj, segname, 'A')
            _renameStrand(strand2Obj, stacks_strand2, ffObj, segname, 'B')

            # connect the two strands
            strandObj = strand1Obj.merge(strand2Obj)
            strandObj._calcZMatrix = False

        # single-strand
        else:
            _renameStrand(strand1Obj, stacks_strand1, ffObj, segname, 'A')

        try:
            strandObj.title = "5'-{}-3'".format(strand1)+' | '+"3'-{}-5'".format(strand2)
            fileutils.archiveFile('double-strand.xyz')
            fileutils.archiveFile('double-strand.pdb')
            strandObj.save('double-strand.xyz')
            strandObj.save('double-strand.pdb')
        except:
            fileutils.archiveFile('single-strand.xyz')
            fileutils.archiveFile('single-strand.pdb')
            strand1Obj.title = "5'-{}-3'".format(strand1)
            strand1Obj.save('single-strand.xyz')
            strand1Obj.save('single-strand.pdb')

        # print a summary
        print("\n >>> Polynucleotide fragment generated")
        print("\n sequence: 5'-{}-3'".format(strand1))
        if strand2:
            print("           3'-{}-5'".format(strand2))
        if backbone:
            print("\n structure reference      : {}".format(backbone))
        print("\n alpha   (C5'-O5'-P-O3')  : {:6.1f} degrees".format(alpha))
        print(" beta    (C4'-C5'-O5'-P)  : {:6.1f} degrees".format(beta))
        print(" gamma   (C3'-C4'-C5'-O5'): {:6.1f} degrees".format(gamma))
        print(" delta   (O3'-C3'-C4'-C5'): {:6.1f} degrees".format(delta))
        print(" epsilon (P-O3'-C3'-C4')  : {:6.1f} degrees".format(epsilon))
        print(" zeta    (O5'-P-O3'-C3')  : {:6.1f} degrees".format(zeta))
        print(" chi     (C2-N1-C1'-O4'")
        print("       or C4-N9-C1'-O4')  : {:6.1f} degrees".format(chi))
        print("\n bases per turn           : {:6.1f}".format(360.0/(rad2deg(arccos(dot(PO0_strand1,PO1_strand1)/linalg.norm(PO0_strand1)/linalg.norm(PO1_strand1))))))
        print(" total length along axis  : {:6.1f} Ångstrom".format(len_strand1*unitconvert.Bohr2Angstrom))
        print(" rise/base along axis     : {:6.1f} Ångstrom".format(len_strand1/nblocks*unitconvert.Bohr2Angstrom))
        print(" P-P diameter of strand   : {:6.1f} Ångstrom\n".format(diameter1*unitconvert.Bohr2Angstrom))

        # save
        filename = 'strand1{}'.format({'a':'_conform_A'}.get(conformation.lower(), ''))

        

        try:
            return strandObj
        except:
            return strand1Obj




class Protein():
    '''Template class for protein'''

    def __new__(cls,
               *args,
                id='',
                source='rcsb',
              **kwargs):
        '''Initialise'''

        from urllib.request import urlretrieve
        from .fragment      import Fragment

        selectcases = { 'rcsb':'https://files.rcsb.org/download/{}.pdb',
                        'ebi' :'https://www.ebi.ac.uk/pdbe/entry-files/download/pdb{}.ent'
                      }

        # TODO
        reports = { 'rcsb':'https://files.rcsb.org/pub/pdb/validation_reports/{}/{}/{}_full_validation.pdf',
                  }

        if id:

            print('\n >>> Downloading protein {} from {}...'.format(id.upper(), source.upper()))

            # PDB url
            url = selectcases.get(source.lower(), 'rcsb').format(id.lower())
            filename = '{}_{}.pdb'.format(id.upper(), source.upper())
            urlretrieve(url, filename)

            inst = Fragment(coords=filename)
            inst.title = id.upper()

        return inst
