TOP.py

##################
## 7 # TOPOLOGY ##  -> @TOP <-
##################
import IO,FUNC,MAP
import logging,math

# This is a generic class for Topology Bonded Type definitions
class Bonded:
    # The init method is generic to the bonded types,
    # but may call the set method if atoms are given
    # as (ID, ResidueName, SecondaryStructure) tuples
    # The set method is specific to the different types.
    def __init__(self,other=None,options=None,**kwargs):
        self.atoms = []
        self.type = -1
        self.parameters = []
        self.comments = []
        self.category = None 

        if options and type(options) == dict:
            self.options = options
        if other:
            # If other is given, then copy the attributes
            # if it is of the same class or set the 
            # attributes according to the key names if
            # it is a dictionary
            if other.__class__ == self.__class__:
                for attr in dir(other):
                    if not attr[0] == "_":
                        setattr(self,attr,getattr(other,attr))
            elif type(other) == dict:
                for attr in other.keys():
                    setattr(self,attr,other[attr])
            elif type(other) in (list,tuple):
                self.atoms = other

        # For every item in the kwargs keys, set the attribute
        # with the same name. This can be used to specify the 
        # attributes directly or to override attributes 
        # copied from the 'other' argument.
        for key in kwargs:
            setattr(self,key,kwargs[key])
        
        # If atoms are given as tuples of
        # (ID, ResidueName[, SecondaryStructure])
        # then determine the corresponding parameters 
        # from the lists above
        if self.atoms and type(self.atoms[0]) == tuple:
            self.set(self.atoms,**kwargs)          

    def __nonzero__(self):
        return bool(self.atoms) and bool(self.parameters) 

    def __str__(self):
        if not self.atoms or not self.parameters:
            return ""
        s = ["%5d" % i for i in self.atoms]
        # For exclusions, no type is defined, which equals -1
        if self.type != -1: s.append(" %5d " % self.type)
        # Print integers and floats in proper format and neglect None terms
        s.extend([FUNC.formatString(i) for i in self.parameters if i != None])
        if self.comments:
            s.append(';')
            if type(self.comments) == str:
                s.append(self.comments)
            else:
                s.extend([str(i) for i in self.comments])
        return " ".join(s)

    def __iadd__(self,num):
        self.atoms = tuple([i+int(num) for i in self.atoms])
        return self

    def __add__(self,num):
        out  = self.__class__(self)
        out += num
        return out

    def __eq__(self,other):
        if type(other) in (list,tuple):
            return self.atoms == other
        else:
            return self.atoms == other.atoms and self.type == other.type and self.parameters == other.parameters

    # This function needs to be overridden for descendents
    def set(self,atoms,**kwargs):
        pass


# The set method of this class will look up parameters for backbone beads
# Side chain bonds ought to be set directly, using the constructor
# providing atom numbers, bond type, and parameters
# Constraints are bonds with kb = None, which can be extracted 
# using the category
class Bond(Bonded):
    def set(self,atoms,**kwargs):
        ids,r,ss,ca     = zip(*atoms)     
        self.atoms      = ids
        self.type       = 1
        self.positionCa = ca
        self.comments   = "%s(%s)-%s(%s)" % (r[0],ss[0],r[1],ss[1])
        # The category can be used to keep bonds sorted
        self.category   = kwargs.get("category")
       
        if not self.parameters:
            self.parameters = self.options['ForceField'].bbGetBond(r,ca,ss)
        # If there are more than two parameters given, user HAS TO define the bond type as the first one.
        # This is to allow use of any gromacs bond type.
        if self.parameters and len(self.parameters) > 2:
            self.type, self.parameters = self.parameters[0], self.parameters[1:]
        # Backbone bonds also can be constraints. We could change the type further on, but this is more general.
        # Even better would be to add a new type: BB-Constraint
        if self.parameters and self.parameters[1] == None:
            self.category = 'Constraint'

    # Overriding __str__ method to suppress printing of bonds with Fc of 0
    def __str__(self):
        if len(self.parameters) > 1 and self.parameters[1] == 0:
            return ""
        return Bonded.__str__(self)


# Similar to the preceding class
class Angle(Bonded):
    def set(self,atoms,**kwargs):
        ids,r,ss,ca     = zip(*atoms)
        self.atoms      = ids
        self.type       = 2
        self.positionCa = ca
        self.comments   = "%s(%s)-%s(%s)-%s(%s)" % (r[0],ss[0],r[1],ss[1],r[2],ss[2])
        self.category   = kwargs.get("category")
        
        self.parameters = self.options['ForceField'].bbGetAngle(r,ca,ss)
        # If there are more than two parameters given, user HAS TO define the angle type as the first one.
        # This is to allow use of any gromacs angle type.
        if self.parameters and len(self.parameters) > 2:
            self.type, self.parameters = self.parameters[0], self.parameters[1:]

# Similar to the preceding class
class Vsite(Bonded):
    def set(self,atoms,**kwargs):
        ids,r,ss,ca     = zip(*atoms)
        self.atoms      = ids
        self.type       = 1
        self.positionCa = ca
        self.comments   = "%s"% (r[0])
        self.category   = kwargs.get("category")
        self.parameters = kwargs.get("parameters") 
        if not self.parameters:
            self.parameters = self.options['ForceField'].bbGetVsite(r,ca,ss)
            # This assumes type one virtual sites, COM of two beads
            if self.parameters and len(self.parameters)>1:
                self.atoms += self.parameters[:-1]
                self.parameters = self.parameters[-1:]

# Similar to the preceding class
class Exclusion(Bonded):
    def set(self,atoms,**kwargs):
        ids,r,ss,ca     = zip(*atoms)
        self.atoms      = ids
        self.positionCa = ca
        self.comments   = "%s"% (r[0])
        self.category   = kwargs.get("category")
        self.parameters = kwargs.get("parameters") 
        if not self.parameters:
            self.parameters = self.options['ForceField'].bbGetExclusion(r,ca,ss)

# Similar to the preceding class
class Pair(Bonded):
    def set(self,atoms,**kwargs):
        ids,r,ss,ca     = zip(*atoms)
        self.atoms      = ids
        self.positionCa = ca
        self.comments   = "%s(%s)-%s(%s)"% (r[0],ss[0],r[1],ss[1])
        self.category   = kwargs.get("category")
        self.type       = kwargs.get("type") or 1 
        if not self.parameters:
            self.parameters = self.options['ForceField'].bbGetPair(r,ca,ss)

# Similar to the preceding class
class Dihedral(Bonded):
    def set(self,atoms,**kwargs):
        ids,r,ss,ca     = zip(*atoms)
        self.atoms      = ids
        self.type       = 1
        self.positionCa = ca
        self.comments   = "%s(%s)-%s(%s)-%s(%s)-%s(%s)" % (r[0],ss[0],r[1],ss[1],r[2],ss[2],r[3],ss[3])
        self.category   = kwargs.get("category")
        if not self.parameters:
            self.parameters = self.options['ForceField'].bbGetDihedral(r,ca,ss)
        # If there are more than two parameters given, user HAS TO define the diheral type as the first one.
        # This is to allow use of any gromacs dihedral type.
        if self.parameters and len(self.parameters)>2:
            self.type,self.parameters = self.parameters[0],self.parameters[1:]

        # XXX Move this to force field files to accommodate polbb and dna dihedrals in a uniform way
        # if ''.join(i for i in ss) == 'FFFF':
        #    # Collagen
        #    self.parameters = self.options['ForceField'].bbDihedDictD['F']
        #elif ''.join(i for i in ss) == 'EEEE' and self.options['ExtendedDihedrals']:
        #    # Use dihedrals
        #    self.parameters = self.options['ForceField'].bbDihedDictD['E']
        #elif set(ss).issubset("H123"):
        #    # Helix
        #    self.parameters = self.options['ForceField'].bbDihedDictD['H']
        #else:
        #    self.parameters = None


# This list allows to retrieve Bonded class items based on the category
# If standard, dictionary type indexing is used, only exact matches are
# returned. Alternatively, partial matching can be achieved by setting
# a second 'True' argument. 
class CategorizedList(list):
    def __getitem__(self,tag): 
        if type(tag) == int:
            # Call the parent class __getitem__
            return list.__getitem__(self,tag)

        if type(tag) == str:
            return [i for i in self if i.category == tag]

        if tag[1]:
            return [i for i in self if tag[0] in i.category]
        else:
            return [i for i in self if i.category == tag[0]]


class Topology:
    def __init__(self,other=None,options=None,name=""):
        self.name        = ''
        self.nrexcl      = 1
        self.atoms       = CategorizedList()
        self.pairs       = CategorizedList()
        self.vsites      = CategorizedList() 
        self.exclusions  = CategorizedList() 
        self.bonds       = CategorizedList()
        self.angles      = CategorizedList()
        self.dihedrals   = CategorizedList()
        self.impropers   = CategorizedList()
        self.constraints = CategorizedList()
        self.posres      = CategorizedList()
        self.sequence    = []
        self.secstruc    = ""
        # Okay, this is sort of funny; we will add a 
        #   #define mapping virtual_sitesn
        # to the topology file, followed by a header
        #   [ mapping ]
        self.mapping     = []
        # For multiscaling we have to keep track of the number of 
        # real atoms that correspond to the beads in the topology
        self.natoms      = 0        
        self.multiscale  = options['multi']

        if options:
            self.options = options
        else:
            self.options = {}

        if not other:
            # Returning an empty instance
            return
        elif isinstance(other,Topology):
            for attrib in ["atoms","vsites","bonds","angles","dihedrals","impropers","constraints","posres"]:
                setattr(self,attrib,getattr(other,attrib,[]))
        elif isinstance(other,IO.Chain):
            if other.type() == "Protein" and other.options['ForceField'].name in ['polbb']:
                self.fromAminoAcidSequencePolBB(other)
            elif other.type() == "Protein":
                self.fromAminoAcidSequence(other)
            elif other.type() == "Nucleic":
                self.fromNucleicAcidSequence(other)
            else:
                # This chain should not be polymeric, but a collection of molecules
                # For each unique residue type fetch the proper moleculetype 
                self.fromMoleculeList(other)
        if name:
            self.name = name

    def __iadd__(self,other):
        if not isinstance(other,Topology):
            other = Topology(other)
        shift     = len(self.atoms)
        last      = self.atoms[-1]
        atoms     = zip(*other.atoms)
        atoms[0]  = [i+shift for i in atoms[0]]   # Update atom numbers
        atoms[2]  = [i+last[2] for i in atoms[2]] # Update residue numbers
        atoms[5]  = [i+last[5] for i in atoms[5]] # Update charge group numbers
        atoms     = zip(*atoms)
        # The zippings above kills all atoms with specified mass (9 long lists) 
        # Let's add some band aid to fix it...
        # This of course doesn't work if there was a secondary structure to keep
        # but think it affects at this stage only the output comment.
        for i in range(len(atoms)):
            if atoms[i][-1] == 0:
                atoms[i] = atoms[i] + ('c',)
        self.atoms.extend(atoms)
        for attrib in ["bonds","vsites","angles","dihedrals","impropers","constraints","posres"]:
            getattr(self,attrib).extend([source+shift for source in getattr(other,attrib)])
        return self

    def __add__(self,other):
        out = Topology(self)
        if not isinstance(other,Topology):
            other = Topology(other)
        out += other
        return out

    def __str__(self):
        if self.multiscale:
             out  = [ '; MARTINI (%s) Multiscale virtual sites topology section for "%s"' %(self.options['ForceField'].name,self.name) ]
        else:
             string  = '; MARTINI (%s) Coarse Grained topology file for "%s"' %(self.options['ForceField'].name, self.name)
             string += '\n; Created by martinize.py version %s \n; Using the following options:  ' %(self.options['Version'])
             string += ' '.join(self.options['Arguments'])
             string += '\n; #####################################################################################################'
             string += '\n; This topology is based on development beta of Martini DNA and should NOT be used for production runs.'
             string += '\n; #####################################################################################################'
             out  = [ string ]
        if self.sequence:
            out += [
                '; Sequence:',
                '; ' + ''.join([ MAP.AA321.get(AA) for AA in self.sequence ]),
                '; Secondary Structure:',
                '; ' + self.secstruc,
                ]
        
        # Do not print a molecule name when multiscaling
        # In that case, the topology created here needs to be appended
        # at the end of an atomistic moleculetype
        if not self.multiscale:
            out += [ '\n[ moleculetype ]',
                     '; Name         Exclusions',  
                     '%-15s %3d' % (self.name,self.nrexcl)]

        out.append('\n[ atoms ]')

        # For virtual sites and dummy beads we have to be able to specify the mass.
        # Thus we need two different format strings:
        fs8 = '%5d %5s %5d %5s %5s %5d %7.4f ; %s'  
        fs9 = '%5d %5s %5d %5s %5s %5d %7.4f %7.4f ; %s'  
        out.extend([len(i)==9 and fs9%i or fs8%i for i in self.atoms])
        #for i in self.atoms:
        #    print i

        # Print the pairs.
        pairs = [str(i) for i in self.pairs]
        if pairs:
            out.append('\n[ pairs ]')
            out.extend(pairs)

        # Print out the vsites only if they excist. Right now it can only be type 1 virual sites.
        # TODO: This needs to be generalized for all virtual site types.
        vsitesBB = [str(i) for i in self.vsites["BB"]]
        vsitesSC = [str(i) for i in self.vsites["SC"]]
        if vsitesBB or vsitesSC:
            out.append('\n[ virtual_sites3 ]')
        if vsitesBB:
            out.append("; Backbone virtual sites.")
            out.extend(vsitesBB)
        if vsitesSC:
            out.append("; Backbone virtual sites.")
            out.extend(vsitesSC)

        # Bonds in order: backbone, backbone-sidechain, sidechain, short elastic, long elastic        
        out.append("\n[ bonds ]")       
        # Backbone-backbone
        bonds = [str(i) for i in self.bonds["BB"]]
        if bonds:
            out.append("; Backbone bonds")
            out.extend(bonds)
        # Rubber Bands
        bonds = [str(i) for i in self.bonds["Rubber",True]]
        #print bonds
        if bonds:
            # Add a CPP style directive to allow control over the elastic network
            out.append("#ifdef RUBBER_BANDS")
            out.append("#ifndef RUBBER_FC\n#define RUBBER_FC %f\n#endif"%self.options['ElasticMaximumForce'])
            out.extend(bonds)
            out.append("#endif")
        # Backbone-Sidechain/Sidechain-Sidechain
        bonds = [str(i) for i in self.bonds["SC"]]
        if bonds:
            out.append("; Sidechain bonds")
            out.extend(bonds)
        # Short elastic/Long elastic
        bonds = [str(i) for i in self.bonds["Elastic short"]]
        if bonds:
            out.append("; Short elastic bonds for extended regions")
            out.extend(bonds)
        bonds = [str(i) for i in self.bonds["Elastic long"]]
        if bonds:
            out.append("; Long elastic bonds for extended regions")
            out.extend(bonds)
        # Cystine bridges
        bonds = [str(i) for i in self.bonds["Cystine"]]
        if bonds:
            out.append("; Cystine bridges")
            out.extend(bonds)
        # Other links
        bonds = [str(i) for i in self.bonds["Link"]]
        if bonds:
            out.append("; Links/Cystine bridges")
            out.extend(bonds)

        # Constraints
        out.append("\n[ constraints ]")
        out.extend([str(i) for i in self.bonds["Constraint"]])

        # Print out the exclusions only if they excist.
        exclusions = [str(i) for i in self.exclusions]
        if exclusions:
            out.append('\n[ exclusions ]')
            out.extend(exclusions)

        if self.multiscale:
            out += ['\n;\n; Coarse grained to atomistic mapping\n;',
                    '#define mapping virtual_sitesn',
                    '[ mapping ]']
            for i,j in self.mapping:
                out.append( ("%5d     2 "%i)+" ".join(["%5d"%k for k in j]) )
            
            logging.info('Created virtual sites section for multiscaled topology')
            return "\n".join(out)

        # Angles
        out.append("\n[ angles ]")
        out.append("; Backbone angles")
        out.extend([str(i) for i in self.angles["BBB"]])
        out.append("; Backbone-sidechain angles")
        out.extend([str(i) for i in self.angles["BBS"]])
        out.append("; Sidechain angles")
        out.extend([str(i) for i in self.angles["SC"]])

        # Dihedrals
        out.append("\n[ dihedrals ]")
        out.append("; Backbone dihedrals")
        out.extend([str(i) for i in self.dihedrals["BBBB"] if i.parameters])
        out.append("; Sidechain dihedrals")
        out.extend([str(i) for i in self.dihedrals["BSC"] if i.parameters])
        out.append("; Sidechain improper dihedrals")
        out.extend([str(i) for i in self.dihedrals["SC"] if i.parameters])

        # Postition Restraints
        if self.posres:
            out.append("\n#ifdef POSRES")
            out.append("#ifndef POSRES_FC\n#define POSRES_FC %.2f\n#endif"%self.options['PosResForce'])
            out.append(" [ position_restraints ]")
            out.extend(['  %5d    1    POSRES_FC    POSRES_FC    POSRES_FC'%i for i in self.posres])
            out.append("#endif")

        logging.info('Created coarsegrained topology')
        return "\n".join(out)


    def ndx(self,start):
        bonds_out = []
        angles_out = []
        dihs_out = []
        # We print out an index file with similar bonds grouped together 
        # We print out an index group for each bond, angle and dihedral
        # Backbone-backbone
        bonds = [str(i) for i in self.bonds["BB"]]
        for i in bonds:
            params = i.split()
            params = [str(start + int(i)) for i in params[:2]]
            bonds_out.append('[BB-bond-%s-%s]\n' % (params[0], params[1]))
            bonds_out.append(' %s %s\n' % (params[0], params[1]))
        # Backbone-Sidechain/Sidechain-Sidechain
        bonds = [str(i) for i in self.bonds["SC"]]
        for i in bonds:
            params = i.split()
            params = [str(start + int(i)) for i in params[:2]]
            bonds_out.append('[SC-bond-%s-%s]\n' % (params[0], params[1]))
            bonds_out.append(' %s %s\n' % (params[0], params[1]))
        # Constraints
        bonds = [str(i) for i in self.bonds["Constraint"]]
        for i in bonds:
            params = i.split()
            params = [str(start + int(i)) for i in params[:2]]
            bonds_out.append('[Const-bond-%s-%s]\n' % (params[0], params[1]))
            bonds_out.append(' %s %s\n' % (params[0], params[1]))

        # Angles
        angles = [str(i) for i in self.angles["BBB"]]
        for i in angles:
            params = i.split()
            params = [str(start + int(i)) for i in params[:3]]
            angles_out.append('[BBB-angle-%s-%s-%s]\n' % (params[0], params[1], params[2]))
            angles_out.append(' %s %s %s\n' % (params[0], params[1], params[2]))
        angles = [str(i) for i in self.angles["BBS"]]
        for i in angles:
            params = i.split()
            params = [str(start + int(i)) for i in params[:3]]
            angles_out.append('[BBS-angle-%s-%s-%s]\n' % (params[0], params[1], params[2]))
            angles_out.append(' %s %s %s\n' % (params[0], params[1], params[2]))
        angles = [str(i) for i in self.angles["SC"]]
        for i in angles:
            params = i.split()
            params = [str(start + int(i)) for i in params[:3]]
            angles_out.append('[SC-angle-%s-%s-%s]\n' % (params[0], params[1], params[2]))
            angles_out.append(' %s %s %s\n' % (params[0], params[1], params[2]))

        # Dihedrals
        dihs = [str(i) for i in self.dihedrals["BBBB"] if i.parameters]
        for i in dihs:
            params = i.split()
            params = [str(start + int(i)) for i in params[:4]]
            dihs_out.append('[BB-dihedral-%s-%s-%s-%s]\n' % (params[0], params[1], params[2], params[3]))
            dihs_out.append(' %s %s %s %s\n' % (params[0], params[1], params[2], params[3]))
        dihs = [str(i) for i in self.dihedrals["BSC"] if i.parameters]
        for i in dihs:
            params = i.split()
            params = [str(start + int(i)) for i in params[:4]]
            dihs_out.append('[BSC-dihedral-%s-%s-%s-%s]\n' % (params[0], params[1], params[2], params[3]))
            dihs_out.append(' %s %s %s %s\n' % (params[0], params[1], params[2], params[3]))
        dihs = [str(i) for i in self.dihedrals["SC"] if i.parameters]
        for i in dihs:
            params = i.split()
            params = [str(start + int(i)) for i in params[:4]]
            dihs_out.append('[SC-dihedral-%s-%s-%s-%s]\n' % (params[0], params[1], params[2], params[3]))
            dihs_out.append(' %s %s %s %s\n' % (params[0], params[1], params[2], params[3]))

        logging.info('Created index files for bonded parameters.')
        return ''.join(bonds_out), ''.join(angles_out), ''.join(dihs_out)

  
    # The sequence function can be used to generate the topology for 
    # a sequence :) either given as sequence or as chain
    def fromAminoAcidSequence(self,sequence,secstruc=None,links=None,breaks=None,
                              mapping=None,rubber=False,multi=False):
        # Shift for the atom numbers of the atomistic part in a chain 
        # that is being multiscaled
        shift = 0
        # First check if we get a sequence or a Chain instance
        if isinstance(sequence, IO.Chain):
            chain         = sequence
            links         = chain.links
            breaks        = chain.breaks
            # If the mapping is not specified, the actual mapping is taken,
            # used to construct the coarse grained system from the atomistic one.
            # The function argument "mapping" could be used to use a default 
            # mapping scheme in stead, like the mapping for the GROMOS96 force field.
            mapping = mapping           or chain.mapping
            multi   = self.options['multi']  or chain.multiscale
            self.secstruc = chain.sstypes or len(chain)*"C"
            self.sequence = chain.sequence
            # If anything hints towards multiscaling, do multiscaling
            self.multiscale = self.multiscale or chain.multiscale or multi
            if self.multiscale:
                shift        = self.natoms
                self.natoms += len(chain.atoms())
        elif not secstruc:
            # If no secondary structure is provided, set all to coil
            chain         = None
            self.secstruc = len(self.sequence)*"C"
        else:
            # If a secondary structure is provided, use that. chain is none.
            chain         = None
            self.secstruc = secstruc

        logging.debug(self.secstruc)
        logging.debug(self.sequence)

        # Fetch the sidechains
        # Pad with empty lists for atoms, bonds, angles 
        # and dihedrals, and take the first four lists out
        # This will avoid errors for residues for which 
        # these are not defined.

        sc = [(self.options['ForceField'].sidechains[res]+5*[[]])[:5] for res in self.sequence]

        # ID of the first atom/residue
        # The atom number and residue number follow from the last 
        # atom c.q. residue id in the list processed in the topology
        # thus far. In the case of multiscaling, the real atoms need 
        # also be accounted for.
        startAtom = self.natoms + 1 
        startResi = self.atoms and self.atoms[-1][2]+1 or 1

        # Backbone bead atom IDs
        bbid = [startAtom]
        for i in zip(*sc)[0]:
            bbid.append(bbid[-1]+len(i)+1)

        # Calpha positions, to get Elnedyn BBB-angles and BB-bond lengths
        # positionCa = [residue[1][4:] for residue in chain.residues]
        # The old method (line above) assumed no hydrogens: Ca would always be
        # the second atom of the residue. Now we look at the name.
        positionCa = []
        for residue in chain.residues:
            for atom in residue:
                if atom[0] == "CA":
                    positionCa.append(atom[4:])

        # Residue numbers for this moleculetype topology
        resid = range(startResi,startResi+len(self.sequence))     
        
        # This contains the information for deriving backbone bead types,
        # bb bond types, bbb/bbs angle types, and bbbb dihedral types and
        # Elnedyn BB-bondlength BBB-angles
        seqss = zip(bbid,self.sequence,self.secstruc,positionCa)

        # Fetch the proper backbone beads          
        bb = [self.options['ForceField'].bbGetBead(res,typ) for num,res,typ,Ca in seqss]

        # If termini need to be charged, change the bead types
        if not self.options['NeutralTermini']:
            bb[0]  ="Qd"
            bb[-1] = "Qa"

        # If breaks need to be charged, change the bead types 
        if self.options['ChargesAtBreaks']:
            for i in breaks:
                bb[i]   = "Qd"
                bb[i-1] = "Qa"

        # For backbone parameters, iterate over fragments, inferred from breaks
        for i,j in zip([0]+breaks,breaks+[-1]):
            # Extract the fragment
            frg = j==-1 and seqss[i:] or seqss[i:j]

            # Iterate over backbone bonds
            self.bonds.extend([Bond(pair,category="BB",options=self.options,) for pair in zip(frg,frg[1:])])

            # Iterate over backbone angles
            # Don't skip the first and last residue in the fragment
            self.angles.extend([Angle(triple,options=self.options,category="BBB") for triple in zip(frg,frg[1:],frg[2:])])

            # Get backbone quadruples
            quadruples = zip(frg,frg[1:],frg[2:],frg[3:])

            # No i-1,i,i+1,i+2 interactions defined for Elnedyn
            if self.options['ForceField'].UseBBBBDihedrals:
                # Process dihedrals
                for q in quadruples:
                    id,rn,ss,ca = zip(*q)
                    # Maybe do local elastic networks
                    if ss == ("E","E","E","E") and not self.options['ExtendedDihedrals']:
                        # This one may already be listed as the 2-4 bond of a previous one
                        if not (id[0],id[2]) in self.bonds:
                            self.bonds.append(Bond(options=self.options,atoms=(id[0],id[2]),parameters=self.options['ForceField'].ebonds['short'],type=1,
                                                   comments="%s(%s)-%s(%s) 1-3"%(rn[0],id[0],rn[2],id[2]),
                                                   category="Elastic short"))
                        self.bonds.append(Bond(options=self.options,atoms=(id[1],id[3]),parameters=self.options['ForceField'].ebonds['short'],type=1,
                                               comments="%s(%s)-%s(%s) 2-4"%(rn[1],id[1],rn[3],id[3]),
                                               category="Elastic short"))
                        self.bonds.append(Bond(options=self.options,atoms=(id[0],id[3]),parameters=self.options['ForceField'].ebonds['long'],type=1,
                                               comments="%s(%s)-%s(%s) 1-4"%(rn[0],id[0],rn[3],id[3]),
                                               category="Elastic long"))
                    else:
                        # Since dihedrals can return None, we first collect them separately and then
                        # add the non-None ones to the list
                        dihed = Dihedral(q,options=self.options,category="BBBB")
                        if dihed:
                            self.dihedrals.append(dihed)

            # Elnedyn does not use backbone-backbone-sidechain-angles
            if self.options['ForceField'].UseBBSAngles:
                # Backbone-Backbone-Sidechain angles
                # If the first residue has a sidechain, we take SBB, otherwise we skip it
                # For other sidechains, we 'just' take BBS
                if len(frg) > 1 and frg[1][0]-frg[0][0] > 1:
                    self.angles.append(Angle(options=self.options,atoms=(frg[0][0]+1,frg[0][0],frg[1][0]),parameters=self.options['ForceField'].bbsangle,type=2,
                                            comments="%s(%s)-%s(%s) SBB"%(frg[0][1],frg[0][2],frg[1][1],frg[1][2]),
                                            category="BBS"))
    
                # Start from first residue: connects sidechain of second residue
                for (ai,ni,si,ci),(aj,nj,sj,cj),s in zip(frg[0:],frg[1:],sc[1:]):
                    if s[0]:
                        self.angles.append(Angle(options=self.options,atoms=(ai,aj,aj+1),parameters=self.options['ForceField'].bbsangle,type=2,
                                                comments="%s(%s)-%s(%s) SBB"%(ni,si,nj,sj),
                                                category="BBS"))
           
        # Now do the atom list, and take the sidechains along
        #
        # AtomID AtomType ResidueID ResidueName AtomName ChargeGroup Charge ; Comments
        # 
        atid  = startAtom
        for resi,resname,bbb,sidechn,ss in zip(resid,self.sequence,bb,sc,self.secstruc):
            scatoms, bon_par, ang_par, dih_par, vsite_par = sidechn

            # Side chain bonded terms
            # Collect bond, angle and dihedral connectivity
            bon_con,ang_con,dih_con,vsite_con = (self.options['ForceField'].connectivity[resname]+4*[[]])[:4]

            # Side Chain Bonds/Constraints
            for atids,par in zip(bon_con,bon_par):
                if par[1] == None:
                    self.bonds.append(Bond(options=self.options,atoms=atids,parameters=[par[0]],type=1,
                                           comments=resname,category="Constraint"))
                else:
                    self.bonds.append(Bond(options=self.options,atoms=atids,parameters=par,type=1,
                                           comments=resname,category="SC"))
                # Shift the atom numbers
                self.bonds[-1] += atid

            # Side Chain Angles
            for atids,par in zip(ang_con,ang_par):
                self.angles.append(Angle(options=self.options,atoms=atids,parameters=par,type=2,
                                         comments=resname,category="SC"))
                # Shift the atom numbers
                self.angles[-1] += atid

            # Side Chain Dihedrals
            for atids,par in zip(dih_con,dih_par):
                self.dihedrals.append(Dihedral(options=self.options,atoms=atids,parameters=par,type=2,
                                               comments=resname,category="SC"))
                # Shift the atom numbers
                self.dihedrals[-1] += atid

            # Side Chain V-Sites
            for atids,par in zip(vsite_con,vsite_par):
                self.vsites.append(Vsite(options=self.options,atoms=atids,parameters=par,type=1,
                                               comments=resname,category="SC"))
                # Shift the atom numbers
                self.vsites[-1] += atid
            
            # Side Chain exclusions
            # The new polarizable forcefield give problems with the charges in the sidechain, if the backbone is also charged.
            # To avoid that, we add explicit exclusions
            if bbb in self.options['ForceField'].charges.keys() and resname in self.options['ForceField'].mass_charge.keys():
                for i in [i for i, d in enumerate(scatoms) if d=='D']:
                    self.exclusions.append(Exclusion(options=self.options,atoms=(atid,i+atid+1),comments='%s(%s)'%(resname,resi),parameters=(None,)))

            # All residue atoms
            counter = 0  # Counts over beads
            for atype,aname in zip([bbb]+list(scatoms),MAP.CoarseGrained.residue_bead_names):
                if self.multiscale:
                    atype,aname = "v"+atype,"v"+aname
                # If mass or charge diverse, we adopt it here. 
                # We don't want to do this for BB beads because of charged termini.
                if resname in self.options['ForceField'].mass_charge.keys() and counter != 0:
                    M,Q = self.options['ForceField'].mass_charge[resname]
                    aname = Q[counter-1]>0 and 'SCP' or Q[counter-1]<0 and 'SCN' or aname
                    self.atoms.append((atid,atype,resi,resname,aname,atid,Q[counter-1],M[counter-1],ss))
                else:
                    self.atoms.append((atid,atype,resi,resname,aname,atid,self.options['ForceField'].charges.get(atype,0),ss))
                # Doing this here save going over all the atoms onesmore.
                # Generate position restraints for all atoms or Backbone beads only.
                if 'all' in self.options['PosRes']:
                    self.posres.append((atid)) 
                elif aname in self.options['PosRes']:
                    self.posres.append((atid))
                if mapping:
                    self.mapping.append((atid,[i+shift for i in mapping[counter]]))
                atid    += 1
                counter += 1

        # The rubber bands are best applied outside of the chain class, as that gives
        # more control when chains need to be merged. The possibility to do it on the 
        # chain level is retained to allow building a complete chain topology in 
        # a straightforward manner after importing this script as module.
        if rubber and chain:
            rubberList = rubberBands(
                [(i[0],j[4:7]) for i,j in zip(self.atoms,chain.cg()) if i[4] in ElasticBeads],
                ElasticLowerBound,ElasticUpperBound,
                ElasticDecayFactor,ElasticDecayPower,
                ElasticMaximumForce,ElasticMinimumForce)
            self.bonds.extend([Bond(i,options=self.options,type=6,category="Rubber band") for i in rubberList])
        
        # Note the equivalent of atomistic atoms that have been processed 
        if chain and self.multiscale:
            self.natoms += len(chain.atoms())

    def fromNucleicAcidSequence(self,sequence,secstruc=None,links=None,breaks=None,
                              mapping=None,rubber=False,multi=False):

        # Shift for the atom numbers of the atomistic part in a chain 
        # that is being multiscaled
        shift = 0
        # Elastic network for DNA
        #self.options['ElasticBeads'] = ['BB1', 'BB2', 'BB3', 'SC1', 'SC2', 'SC3', 'SC4']
        #self.options['ElasticBeads'] = ['BB3', 'SC1', 'SC2', 'SC3', 'SC4']
        #self.options['ElasticBeads'] = ['BB1', 'BB2', 'BB3']
        #self.options['ElasticBeads'] = ['BB1', 'BB2']
        # First check if we get a sequence or a Chain instance
        if isinstance(sequence, IO.Chain):
            chain         = sequence
            links         = chain.links
            breaks        = chain.breaks
            # If the mapping is not specified, the actual mapping is taken,
            # used to construct the coarse grained system from the atomistic one.
            # The function argument "mapping" could be used to use a default 
            # mapping scheme in stead, like the mapping for the GROMOS96 force field.
            mapping = mapping           or chain.mapping
            multi   = self.options['multi']  or chain.multiscale
            self.secstruc = chain.sstypes or len(chain)*"C"
            self.sequence = chain.sequence
            # If anything hints towards multiscaling, do multiscaling
            self.multiscale = self.multiscale or chain.multiscale or multi
            if self.multiscale:
                shift        = self.natoms
                self.natoms += len(chain.atoms())
        elif not secstruc:
            # If no secondary structure is provided, set all to coil
            chain         = None
            self.secstruc = len(self.sequence)*"C"
        else:
            # If a secondary structure is provided, use that. chain is none.
            chain         = None
            self.secstruc = secstruc

        logging.debug(self.secstruc)
        logging.debug(self.sequence)

        # Fetch the base information 
        # Pad with empty lists for atoms, bonds, angles 
        # and dihedrals, and take the first five lists out
        # This will avoid errors for residues for which 
        # these are not defined.

        sc = [(self.options['ForceField'].bases[res]+8*[[]])[:8] for res in self.sequence]

        # ID of the first atom/residue
        # The atom number and residue number follow from the last 
        # atom c.q. residue id in the list processed in the topology
        # thus far. In the case of multiscaling, the real atoms need 
        # also be accounted for.
        startAtom = self.natoms + 1 
        startResi = self.atoms and self.atoms[-1][2]+1 or 1

        # Backbone bead atom IDs
        bbid = [[startAtom,startAtom+1,startAtom+2]]
        for i in zip(*sc)[0]:
            bbid1 = bbid[-1][0]+len(i)+3
            bbid.append([bbid1,bbid1+1,bbid1+2])
            #bbid.append(bbid[-1]+len(i)+1)

        # Residue numbers for this moleculetype topology
        resid = range(startResi,startResi+len(self.sequence))     

        # This contains the information for deriving backbone bead types,
        # bb bond types, bbb/bbs angle types, and bbbb dihedral types.
        seqss = zip(bbid,self.sequence,self.secstruc)

        # Fetch the proper backbone beads          
        # Since there are three beads we need to split these to the list
        bb = [self.options['ForceField'].bbGetBead(res,typ) for num,res,typ in seqss]
        bbMulti = [i for j in bb for i in j]

        # This is going to be usefull for the type of the last backbone bead.
        # If termini need to be charged, change the bead types
        #if not self.options['NeutralTermini']:
        #    bb[0]  ="Qd"
        #    bb[-1] = "Qa"

        # If breaks need to be charged, change the bead types 
        #if self.options['ChargesAtBreaks']:
        #    for i in breaks:
        #        bb[i]   = "Qd"
        #        bb[i-1] = "Qa"

        # For backbone parameters, iterate over fragments, inferred from breaks
        for i,j in zip([0]+breaks,breaks+[-1]):
            # Extract the fragment
            frg = j==-1 and seqss[i:] or seqss[i:j]
            # Expand the 3 bb beads per residue into one long list
            # Resulting list contains three tuples per residue 
            # We use the useless ca parameter to get the correct backbone bond from bbGetBond 
            # The i is used to define ca because it conveniently here gives the right sequence 0,1,2,0,1,...
            frg = [(j[0][i],j[1],j[2],i) for j in frg for i in range(len(j[0]))]

            # Iterate over backbone bonds, two loops are needed because there are multipe cross bonds in the BB.
            # Since bonded interactions can return None, we first collect them separately and then
            # add the non-None ones to the list
            # This part assumes you have a linear backbone (checks only combinations that make sense for such) 
            for ind,k in enumerate(frg):
                    
                # Exclusions iterate over the second atom.
                # and then checked for each pair.
                for l in frg[ind:ind+3]:
                    excl = Exclusion((k,l),category="BB",options=self.options,)
                    if excl:
                        self.exclusions.append(excl) 

                # Pairs iterate over the second atom.
                # and then checked for each pair.
                for l in frg[ind:ind+3]:
                    pair = Pair((k,l),category="BB",options=self.options,)
                    if pair:
                        self.pairs.append(pair) 

                # Bonds iterate over the second atom.
                for l in frg[ind:ind+3]:
                    bond = Bond((k,l),category="BB",options=self.options,)
                    if bond:
                        self.bonds.append(bond)
                        
                    # The angles need a third loop.
                    for m in frg[ind:ind+3]:
                        angle = Angle((k,l,m,),options=self.options,category="BBB")
                        if angle:
                            self.angles.append(angle)
                            
                        # The dihedrals need a fourth loop.
                        for n in frg[ind+1:ind+4]:
                            dihed = Dihedral((k,l,m,n,),options=self.options,category="BBBB")
                            if dihed:
                                if len(dihed.parameters) > 3:
                                    dihed2 = Dihedral((k,l,m,n,),options=self.options, parameters=dihed.parameters[3:],category="BBBB")
                                    dihed.parameters = dihed.parameters[:3]
                                    self.dihedrals.append(dihed)
                                    self.dihedrals.append(dihed2)
                                else:
                                    self.dihedrals.append(dihed)

        # Now do the atom list, and take the sidechains along
        #
        atid  = startAtom
        # We need to do some trickery to get all 3 bb beads in to these lists
        # This adds each element to a list three times, feel free to shorten up
        residMulti = [i for i in resid for j in range(3)]
        sequenceMulti = [i for i in self.sequence for j in range(3)]
        scMulti = [i for i in sc for j in range(3)]
        secstrucMulti = [i for i in self.secstruc for j in range(3)]
        count = 0
        for resi,resname,bbb,sidechn,ss in zip(residMulti,sequenceMulti,bbMulti,scMulti,secstrucMulti):
            # We only want one side chain per three backbone beads so this skips the others
            if (count % 3) == 0:    
                # Note added impropers in contrast to aa
                scatoms, bon_par, ang_par, dih_par, imp_par, vsite_par, excl_par, pair_par = sidechn

                # Side chain bonded terms
                # Collect bond, angle and dihedral connectivity
                # Impropers needed to be added here for DNA
                bon_con,ang_con,dih_con,imp_con,vsite_con,excl_con, pair_con = (self.options['ForceField'].base_connectivity[resname]+7*[[]])[:7]

                # Side Chain Bonds/Constraints
                for atids,par in zip(bon_con,bon_par):
                    if par[2] == None:
                        self.bonds.append(Bond(options=self.options,atoms=atids,parameters=[par[1]],type=par[0],
                                               comments=resname,category="Constraint"))
                    else:
                        self.bonds.append(Bond(options=self.options,atoms=atids,parameters=par[1:],type=par[0],
                                               comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.bonds[-1] += atid

                # Side Chain Angles
                for atids,par in zip(ang_con,ang_par):
                    self.angles.append(Angle(options=self.options,atoms=atids,parameters=par[1:],type=par[0],
                                             comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.angles[-1] += atid

                # Side Chain Dihedrals
                for atids,par in zip(dih_con,dih_par):
                    self.dihedrals.append(Dihedral(options=self.options,atoms=atids,parameters=par[1:],type=par[0],
                                                   comments=resname,category="BSC"))
                    # Shift the atom numbers
                    self.dihedrals[-1] += atid

                # Side Chain Impropers
                for atids,par in zip(imp_con,imp_par):
                    self.dihedrals.append(Dihedral(options=self.options,atoms=atids,parameters=par[1:],type=par[0],
                                                   comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.dihedrals[-1] += atid

                # Side Chain V-Sites
                for atids,par in zip(vsite_con,vsite_par):
                    self.vsites.append(Vsite(options=self.options,atoms=atids,parameters=par,type=1,
                                                   comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.vsites[-1] += atid

                # XXX Should not be hard coded here!
                # Currently DNA needs exclusions for the base
                # The loop runs over the first backbone bead so 3 needs to be added to the indices
                # for i in range(len(scatoms)):
                #    for j in range(i+1, len(scatoms)):
                #        self.exclusions.append(Exclusion(options=self.options,atoms=(i+atid+3,j+atid+3),comments='%s(%s)'%(resname,resi),parameters=(None,)))
                # The above one is ugly hack if only one connectivity record can be used
                # The one below is equally ugly but works for now
                for atids,par in zip(excl_con,excl_par):
                    self.exclusions.append(Exclusion(options=self.options,atoms=atids,parameters=' '))
                    # Shift the atom numbers
                    self.exclusions[-1] += atid

                # Pairs
                for atids,par in zip(pair_con,pair_par):
                    self.pairs.append(Pair(options=self.options,atoms=atids,parameters=' '))
                    # Shift the atom numbers
                    self.pairs[-1] += atid
                
                # All residue atoms
                counter = 0  # Counts over beads
                # Need to tweak this to get all the backbone beads to the list with the side chain
                bbbset = [bbMulti[count], bbMulti[count+1], bbMulti[count+2]]
                for atype,aname in zip(bbbset+list(scatoms),MAP.CoarseGrained.residue_bead_names_dna):
                    if self.multiscale:
                        atype,aname = "v"+atype,"v"+aname
                    #self.atoms.append((atid,atype,resi,resname,aname,atid,self.options['ForceField'].getCharge(atype,aname),ss))
                    if atid in [vSite.atoms[0] for vSite in self.vsites]:
                        #print 'in vsites', atid
                        charge = self.options['ForceField'].getCharge(atype,aname)
                        mass = 0
                        self.atoms.append((atid,atype,resi,resname,aname,atid,charge,mass,ss))
                        #print self.atoms[-1]
                    else:
                        charge = self.options['ForceField'].getCharge(atype,aname)
                        if aname == 'BB1':
                            mass = 72
                        # For virtual sites
                        #elif 'SC' in aname and (resname == 'DA' or resname == 'DG'):
                        #    mass = 60
                        else: 
                            mass = 45
                        self.atoms.append((atid,atype,resi,resname,aname,atid,charge,mass,ss))
                    # Doing this here saves going over all the atoms onesmore.
                    # Generate position restraints for all atoms or Backbone beads only.
                    if 'all' in self.options['PosRes']:
                        self.posres.append((atid)) 
                    elif aname in self.options['PosRes']:
                        self.posres.append((atid))
                    if mapping:
                        self.mapping.append((atid,[i+shift for i in mapping[counter]]))
                    atid    += 1
                    counter += 1
            count += 1

        # One more thing, we need to remove dihedrals (2) and an angle (1)  that reach beyond the 3' end
        # This is stupid to do now but the total number of atoms seems not to be available before
        # This iterate the list in reverse order so that removals don't affect later checks
        for i in range(len(self.dihedrals)-1,-1,-1):
            if (max(self.dihedrals[i].atoms) > self.atoms[-1][0]):
                del self.dihedrals[i]
        for i in range(len(self.angles)-1,-1,-1):
            if (max(self.angles[i].atoms) > self.atoms[-1][0]):
                del self.angles[i]

        # The first residue of DNA (The 5' end) has a backbone that is only two beads. This is ugly
        # but we'll save our headache of changing the whole internal logic. 
        for i in range(len(self.bonds)-1,-1,-1):
            if 1 in self.bonds[i].atoms:
                del self.bonds[i]
            else: 
                self.bonds[i].atoms = tuple([j-1 for j in self.bonds[i].atoms])
        for i in range(len(self.angles)-1,-1,-1):
            if 1 in self.angles[i].atoms:
                del self.angles[i]
            else: 
                self.angles[i].atoms = tuple([j-1 for j in self.angles[i].atoms])
        for i in range(len(self.dihedrals)-1,-1,-1):
            if 1 in self.dihedrals[i].atoms:
                del self.dihedrals[i]
            else: 
                self.dihedrals[i].atoms = tuple([j-1 for j in self.dihedrals[i].atoms])
        for i in range(len(self.vsites)-1,-1,-1):
            if 1 in self.vsites[i].atoms:
                del self.vsites[i]
            else: 
                self.vsites[i].atoms = tuple([j-1 for j in self.vsites[i].atoms])
        for i in range(len(self.exclusions)-1,-1,-1):
            if 1 in self.exclusions[i].atoms:
                del self.exclusions[i]
            else: 
                self.exclusions[i].atoms = tuple([j-1 for j in self.exclusions[i].atoms])
        for i in range(len(self.pairs)-1,-1,-1):
            if 1 in self.pairs[i].atoms:
                del self.pairs[i]
            else: 
                self.pairs[i].atoms = tuple([j-1 for j in self.pairs[i].atoms])

        # Remove the first bead and shift the numbers of all others by one.
        del self.atoms[0]
        for i in range(len(self.atoms)):
            self.atoms[i] = tuple([self.atoms[i][0]-1] + [j for j in self.atoms[i][1:]])
        enStrandLengths.append(self.atoms[-1][0])

        # Remove the last posres, faster than removing first and shifting all.
        if len(self.posres) > 0:
            del self.posres[-1]

        # The rubber bands are best applied outside of the chain class, as that gives
        # more control when chains need to be merged. The possibility to do it on the 
        # chain level is retained to allow building a complete chain topology in 
        # a straightforward manner after importing this script as module.
        if rubber and chain:
            rubberList = rubberBands(
                [(i[0],j[4:7]) for i,j in zip(self.atoms,chain.cg()) if i[4] in ElasticBeads],
                ElasticLowerBound,ElasticUpperBound,
                ElasticDecayFactor,ElasticDecayPower,
                ElasticMaximumForce,ElasticMinimumForce)
            #self.bonds.extend([Bond(i,options=self.options,type=6,category="Rubber band") for i in rubberList])


    def fromAminoAcidSequencePolBB(self,sequence,secstruc=None,links=None,breaks=None,
                              mapping=None,rubber=False,multi=False):

        # Shift for the atom numbers of the atomistic part in a chain 
        # that is being multiscaled
        shift = 0
        # First check if we get a sequence or a Chain instance
        if isinstance(sequence, IO.Chain):
            chain         = sequence
            links         = chain.links
            breaks        = chain.breaks
            # If the mapping is not specified, the actual mapping is taken,
            # used to construct the coarse grained system from the atomistic one.
            # The function argument "mapping" could be used to use a default 
            # mapping scheme in stead, like the mapping for the GROMOS96 force field.
            mapping = mapping           or chain.mapping
            multi   = self.options['multi']  or chain.multiscale
            self.secstruc = chain.sstypes or len(chain)*"C"
            self.sequence = chain.sequence
            # If anything hints towards multiscaling, do multiscaling
            self.multiscale = self.multiscale or chain.multiscale or multi
            if self.multiscale:
                shift        = self.natoms
                self.natoms += len(chain.atoms())
        elif not secstruc:
            # If no secondary structure is provided, set all to coil
            chain         = None
            self.secstruc = len(self.sequence)*"C"
        else:
            # If a secondary structure is provided, use that. chain is none.
            chain         = None
            self.secstruc = secstruc

        logging.debug(self.secstruc)
        logging.debug(self.sequence)

        # Fetch the base information 
        # Pad with empty lists for atoms, bonds, angles 
        # and dihedrals, and take the first five lists out
        # This will avoid errors for residues for which 
        # these are not defined.

        sc = [(self.options['ForceField'].sidechains[res]+6*[[]])[:6] for res in self.sequence]

        # ID of the first atom/residue
        # The atom number and residue number follow from the last 
        # atom c.q. residue id in the list processed in the topology
        # thus far. In the case of multiscaling, the real atoms need 
        # also be accounted for.
        startAtom = self.natoms + 1 
        startResi = self.atoms and self.atoms[-1][2]+1 or 1

        # Backbone bead atom IDs
        # XXX Number of backbone beads hardcoded
        bbid = [[startAtom+i for i in range(4)]]
        for i in zip(*sc)[0]:
            bbid1 = bbid[-1][0]+len(i)+4
            bbid.append([bbid1+i for i in range(4)])
            #bbid.append(bbid[-1]+len(i)+1)

        # Residue numbers for this moleculetype topology
        resid = range(startResi,startResi+len(self.sequence))     

        # This contains the information for deriving backbone bead types,
        # bb bond types, bbb/bbs angle types, and bbbb dihedral types.
        seqss = zip(bbid,self.sequence,self.secstruc)
        # The last residue only has a Calpha in the BB and no dipole.
        seqss[-1] = (seqss[-1][0][:1],)+seqss[-1][1:]

        # Fetch the proper backbone beads          
        # Since there are  N beads we need to split these to the list
        bb = [self.options['ForceField'].bbGetBead(res,typ) for num,res,typ in seqss]
        # bbGetBead always returns all the BB-beads. For the last residue to get out
        # only the last we preprocess it together with seqss (which does have the right length).
        bbMulti = [beadname for residue,beads in zip(seqss,bb) for atomid,beadname in zip(residue[0],beads)]

        # This is going to be usefull for the type of the last backbone bead.
        # If termini need to be charged, change the bead types
        #if not self.options['NeutralTermini']:
        #    bb[0]  ="Qd"
        #    bb[-1] = "Qa"

        # If breaks need to be charged, change the bead types 
        #if self.options['ChargesAtBreaks']:
        #    for i in breaks:
        #        bb[i]   = "Qd"
        #        bb[i-1] = "Qa"

        # For backbone parameters, iterate over fragments, inferred from breaks
        for i,j in zip([0]+breaks,breaks+[-1]):
            # Extract the fragment
            frg = j==-1 and seqss[i:] or seqss[i:j]
            # Expand the N bb beads per residue into one long list
            # Resulting list contains three tuples per residue 
            # We use the useless ca parameter to get the correct backbone bond from bbGetBond 
            frg = [(j[0][i],j[1],j[2],j[0][i]-1) for j in frg for i in range(len(j[0]))]

            # Iterate over backbone bonds, two loops are needed because there are multipe cross bonds in the BB.
            # Since bonded interactions can return None, we first collect them separately and then
            # add the non-None ones to the list
            # Number of backbone beads hardcoded to 4
            for ind,k in enumerate(frg):
                # Virtual site and exclusions we check per bead
                vsite = Vsite((k,),category="BB",options=self.options,)
                if vsite:
                     self.vsites.append(vsite)
                excl = Exclusion((k,),category="BB",options=self.options,)
                if excl and max(excl.parameters) <= len(bbMulti):
                    self.exclusions.append(excl)
                
                # Bonds and pairs interate over the second atom.
                for l in frg[ind:ind+9]:
                    bond = Bond((k,l),category="BB",options=self.options,)
                    if bond:
                        self.bonds.append(bond)
                    pair = Pair((k,l),category="BB",options=self.options,type=2,)
                    if pair:
                        self.pairs.append(pair)

                    # The angles need a third loop.
                    for m in frg[ind:ind+9]:
                        angle = Angle((k,l,m,),options=self.options,category="BBB")
                        if angle:
                            self.angles.append(angle)

                        # The dihedrals need a fourth loop.
                        for n in frg[ind:ind+9]:
                            dihed = Dihedral((k,l,m,n,),options=self.options,category="BBBB")
                            if dihed:
                                self.dihedrals.append(dihed)

        # Now do the atom list, and take the sidechains along
        #
        atid  = startAtom
        # We need to do some trickery to get all N bb beads in to these lists
        # This adds each element to a list three times, feel free to shorten up
        residMulti = [res for i,res in enumerate(resid) for j in range(len(bb[i]))]
        sequenceMulti = [resname for i,resname in enumerate(self.sequence) for j in range(len(bb[i]))]
        # Would it be easier to only give a sidechain for the beads that connect to one (ie the first)?
        scMulti = [sidechain for i,sidechain in enumerate(sc) for j in range(len(bb[i]))]
        secstrucMulti = [secstruc for i,secstruc in enumerate(self.secstruc) for j in range(len(bb[i]))]
        count = 0
        for resi,resname,bbb,sidechn,ss in zip(residMulti,sequenceMulti,bbMulti,scMulti,secstrucMulti):
            # We only want one side chain per three backbone beads so this skips the others
            if (count % len(bb[0])) == 0:    
                # Note added impropers in contrast to aa
                scatoms, bon_par, ang_par, dih_par, imp_par, vsite_par = sidechn

                # Side chain bonded terms
                # Collect bond, angle and dihedral connectivity
                bon_con,ang_con,dih_con,imp_con,vsite_con = (self.options['ForceField'].connectivity[resname]+5*[[]])[:5]

                # Side Chain Bonds/Constraints
                for atids,par in zip(bon_con,bon_par):
                    if par[1] == None:
                        self.bonds.append(Bond(options=self.options,atoms=atids,parameters=[par[0]],type=1,
                                               comments=resname,category="Constraint"))
                    else:
                        self.bonds.append(Bond(options=self.options,atoms=atids,parameters=par,type=1,
                                               comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.bonds[-1] += atid

                # Side Chain Angles
                for atids,par in zip(ang_con,ang_par):
                    self.angles.append(Angle(options=self.options,atoms=atids,parameters=par,type=2,
                                             comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.angles[-1] += atid

                # Side Chain Dihedrals
                for atids,par in zip(dih_con,dih_par):
                    self.dihedrals.append(Dihedral(options=self.options,atoms=atids,parameters=par,type=2,
                                                   comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.dihedrals[-1] += atid

                # Side Chain Impropers
                for atids,par in zip(imp_con,imp_par):
                    self.dihedrals.append(Dihedral(options=self.options,atoms=atids,parameters=par,type=2,
                                                   comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.dihedrals[-1] += atid

                # Side Chain V-Sites
                for atids,par in zip(vsite_con,vsite_par):
                    self.vsites.append(Vsite(options=self.options,atoms=atids,parameters=par,type=1,
                                                   comments=resname,category="SC"))
                    # Shift the atom numbers
                    self.vsites[-1] += atid
            
                # Side Chain exclusions
                # The new polarizable forcefield give problems with the charges in the sidechain, if the backbone is also charged.
                # To avoid that, we add explicit exclusions
                if bbb in self.options['ForceField'].charges.keys() and resname in self.options['ForceField'].mass_charge.keys():
                    for i in [i for i, d in enumerate(scatoms) if d=='D']:
                        self.exclusions.append(Exclusion(options=self.options,atoms=(atid,i+atid+1),comments='%s(%s)'%(resname,resi),parameters=(None,)))

                # All residue atoms
                counter = 0  # Counts over beads
                # Count it the current backbone bead. bbb is the current backbone residue
                bbbset = [bbMulti[count+i] for i in range(len(seqss[resi-1][0]))]
                aNames = MAP.CoarseGrained.residue_bead_names_polBB[1:4]+MAP.CoarseGrained.residue_bead_names_polBB[:1]+MAP.CoarseGrained.residue_bead_names_polBB[4:]
                for atype,aname in zip(bbbset+list(scatoms),aNames):
                    if self.multiscale:
                        atype,aname = "v"+atype,"v"+aname
                    # If mass or charge diverse, we adopt it here. 
                    # We don't want to do this for BB beads because of charged termini.
                    if resname in self.options['ForceField'].mass_charge.keys() and counter != 0:
                        logging.error('Polarizable backbones and sidechains are not yet compatible.')
                        sys.exit()
                        M,Q = self.options['ForceField'].mass_charge[resname]
                        aname = Q[counter-1]>0 and 'SCP' or Q[counter-1]<0 and 'SCN' or aname
                        self.atoms.append((atid,atype,resi,resname,aname,atid,Q[counter-1],M[counter-1],ss))
                    elif atid in [vSite.atoms[0] for vSite in self.vsites]:
                        charge = self.options['ForceField'].getCharge(atype,aname)
                        mass = 0
                        self.atoms.append((atid,atype,resi,resname,aname,atid,charge,mass,ss))
                    else:
                        charge = self.options['ForceField'].getCharge(atype,aname)
                        self.atoms.append((atid,atype,resi,resname,aname,atid,charge,ss))
                    # Doing this here save going over all the atoms onesmore.
                    # Generate position restraints for all atoms or Backbone beads only.
                    if 'all' in self.options['PosRes']:
                        self.posres.append((atid)) 
                    elif aname in self.options['PosRes']:
                        self.posres.append((atid))
                    if mapping:
                        self.mapping.append((atid,[i+shift for i in mapping[counter]]))
                    atid    += 1
                    counter += 1
            count += 1

#        # One more thing, we need to remove dihedrals (2) and an angle (1)  that reach beyond the 3' end
#        # This is stupid to do now but the total number of atoms seems not to be available before
#        # This iterate the list in reverse order so that removals don't affect later checks
#        for i in range(len(self.dihedrals)-1,-1,-1):
#            if (max(self.dihedrals[i].atoms) > self.atoms[-1][0]):
#                del self.dihedrals[i]
#        for i in range(len(self.angles)-1,-1,-1):
#            if (max(self.angles[i].atoms) > self.atoms[-1][0]):
#                del self.angles[i]


    def fromMoleculeList(self,other):
        pass