Skip to content

polymer

Module documentation for polymer.

API Reference

ArbdModel

Bases: PdbModel

Source code in mrdna/arbdmodel/submodule/__init__.py
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
class ArbdModel(PdbModel):

    class _GroupSite():
        """ Class to represent a collection of particles that can be used by bond potentials """
        def __init__(self, particles, weights=None):
            if weights is not None:
                raise NotImplementedError
            self.particles = particles
            self.idx = None
            self.restraints = []

        def get_center(self):
            c = np.array((0,0,0))
            for p in self.particles:
                c = c + p.get_collapsed_position()
            c = c / len(self.particles)
            return c

        def add_restraint(self, restraint):
            self.restraints.append( restraint )
        def get_restraints(self):
            return [(self,r) for r in self.restraints]


    def __init__(self, children, origin=None, dimensions=(1000,1000,1000),
                 remove_duplicate_bonded_terms=True,
                 configuration=None, dummy_types=tuple(), **conf_params):

        PdbModel.__init__(self, children, dimensions, remove_duplicate_bonded_terms)
        self.origin = origin

        if configuration is None: 
            configuration = SimConf(**conf_params)
        self.configuration = configuration

        self.num_particles = 0
        self.particles = []
        self.type_counts = None

        self.dummy_types = dummy_types # TODO, determine whether these are really needed
        if len(self.dummy_types) != 0:
            raise("Dummy types have been deprecated")

        self.nonbonded_interactions = []
        self._nonbonded_interaction_files = [] # This could be made more robust

        self.cacheUpToDate = False

        self.group_sites = []

    def add(self, obj):
        self._clear_types()     # TODO: call upon (potentially nested) `child.add(add)`
        return Parent.add(self, obj)

    def add_group_site(self, particles, weights=None):
        g = ArbdModel._GroupSite(particles, weights)
        self.group_sites.append(g)
        return g

    def _clear_types(self):
        devlogger.debug(f'{self}: Clearing types') 
        self.type_counts = None


    def clear_all(self, keep_children=False):
        Parent.clear_all(self, keep_children=keep_children)
        self.particles = []
        self.num_particles = 0
        self._clear_types()
        self.group_sites = []
        self._nonbonded_interaction_files = []

    def extend(self, other_model, copy=False):

        if any( [p.rigid for p in self] + [p.rigid for p in other_model] ):
            raise NotImplementedError('Models containing rigid bodies cannot yet be extended')

        assert( isinstance(other_model, ArbdModel) )
        if copy == True:
            logger.warning(f'Forcing {self.__class__}.extend(other_model,copy=False)')
            copy = False

        ## Combine particle types, taking care to handle name clashes
        self._countParticleTypes()
        other_model._countParticleTypes()

        self.getParticleTypesAndCounts()
        other_model.getParticleTypesAndCounts()

        names = set([t.name for t in self.type_counts.keys()])

        devlogger.debug(f'Combining types {self.type_counts.keys()} and {other_model.type_counts.keys()}')
        for t1 in self.type_counts.keys():
            for t2 in other_model.type_counts.keys():
                if t1.name == t2.name and t1.__eq__(t2, check_equal=False):
                    i = 1
                    new_name = f'{t2.name}{i}'
                    while new_name in names:
                        i += 1
                        new_name = f'{t2.name}{i}'
                    devlogger.debug(f'Updating {t2.name} to {new_name}')
                    t2.name = new_name
                    names.add(new_name)

        ## Combine interactions
        for i, tA, tB in other_model.nonbonded_interactions:
            devlogger.debug(f'Combining model interactions {i} {tA} {tB}')
            self.add_nonbonded_interaction(i,tA,tB)

        # for g in other_model.children:
        #     self.update(g, copy=copy)
        g = Group()
        for attr in 'children position orientation bonds angles dihedrals impropers exclusions vector_angles bond_angles product_potentials group_sites'.split():
            g.__setattr__(attr, other_model.__getattribute__(attr))
        devlogger.debug(f'Updating {self} with {g.children[0].children[0]}')
        self.update(g, copy=copy)

        self._clear_types()

        ## Combine configurations
        self.configuration = other_model.configuration.combine(self.configuration, policy='best', warn=True)

    def assign_IBI_degrees_of_freedom(self):

        """ Convenience routine that adds degrees of freedom to
        corresponding IBI potentials """

        from .ibi import BondDof, AngleDof, DihedralDof, PairDistributionDof

        self.bonded_ibi_potentials = set()
        self.nonbonded_ibi_potentials = set()

        for get_fn, parts_pot_fn, cls in (
                (self.get_bonds,     lambda x: (x[:2],x[2]), BondDof),
                (self.get_angles,    lambda x: (x[:3],x[3]), AngleDof),
                (self.get_dihedrals, lambda x: (x[:4],x[4]), DihedralDof)
        ):
            for x in get_fn():
                parts, pot = parts_pot_fn(x)
                try: pot.degrees_of_freedom
                except: continue
                pot.degrees_of_freedom.append( cls(*parts) )
                if pot not in self.bonded_ibi_potentials:
                    self.bonded_ibi_potentials.add( pot )

        logger.info(f'Gathering nonbonded IBI degrees of freedom')
        all_exclusions = self.get_exclusions()
        try:
            raise               # particle.idx wasn't set without this
            if len(self.type_counts) == 0: raise
        except:
            self.getParticleTypesAndCounts()

        type_to_particles = dict()
        for p in self:
            t = p.type_
            if t not in type_to_particles: type_to_particles[t] = []
            type_to_particles[t].append(p)
        type_to_particles[None] = [p for p in self]

        already_handled = set()
        for pot, tA, tB in self.nonbonded_interactions:
            ## Avoid including particle type pairs that don't apply
            if (tA,tB) in already_handled: continue
            already_handled.add((tA,tB))
            already_handled.add((tB,tA))

            try:
                pot.degrees_of_freedom
            except:
                continue

            def _cond(pair):
                b1 = (tA is None or pair[0].type_ == tA) and (tB is None or pair[1].type_ == tB)
                b2 = (tB is None or pair[0].type_ == tB) and (tA is None or pair[1].type_ == tA)
                return (b1 or b2)

            ex = list(filter(_cond, all_exclusions))
            dof = PairDistributionDof( type_to_particles[tA], type_to_particles[tB], exclusions=ex,
                                       range_=pot.range_, resolution=pot.resolution)

            pot.degrees_of_freedom.append( dof )

            if pot not in self.nonbonded_ibi_potentials:
                self.nonbonded_ibi_potentials.add( pot )

    def load_target_IBI_distributions(self):
        raise NotImplementedError

    def run_IBI(self, iterations, directory = './', engine = None, replicas = 1, run_minimization = True, first_iteration=1, target_universe = None, target_trajectory = None):
        try:
            assert( len(self.bonded_ibi_potentials) + len(self.nonbonded_ibi_potentials) > 0 )
        except:
            raise ValueError('Model does not appear to contain IBI potentials; perhaps you forgot to run self.assign_IBI_degrees_of_freedom()')
        logger.info(f'Running {iterations} IBI iterations with {len(self.bonded_ibi_potentials)} bonded and {len(self.nonbonded_ibi_potentials)} nonbonded IBI potentials')

        if engine is None:
            engine = ArbdEngine(
                num_steps = 5e6,
                output_period = 1e4,
            )

        if np.array(self.dimensions).size > 3:
            raise NotImplementedError('IBI only implemented for systems with orthorhombic unit cells')
        else:
            box = tuple(list(self.dimensions[:3]) + ([90]*3))

        restart_file = None

        if target_universe is not None:
            _potentials = list(self.bonded_ibi_potentials)+list(self.nonbonded_ibi_potentials)
            with logging_redirect_tqdm(loggers=[logger,devlogger]):
                for potential in tqdm(_potentials, desc='Calculating target distributions'):
                    potential.get_target_distribution(target_universe, trajectory=target_trajectory)

        def _load_cg_u(iteration):
            name = f'ibi-{iteration:03d}'
            psf = '{}/{}.psf'.format(directory,name)
            globstring=f'{directory}/output/{name}.*dcd'
            dcds = [f for f in glob(globstring) if 'momentum' not in f]
            if len(dcds) == 0: raise ValueError(f'Expected to find dcds at {globstring}')
            cg_u = mda.Universe(psf,*dcds)
            return cg_u

        cg_u = None
        if first_iteration > 1:
            cg_u = _load_cg_u(first_iteration-1)
            with logging_redirect_tqdm(loggers=[logger,devlogger]):
                for p in tqdm(_potentials, desc='Calculating initial CG distributions'):
                    p.get_cg_distribution(cg_u, box=box, recalculate=False)

        for p in _potentials:
            p.iteration = first_iteration

        for i in range(first_iteration, iterations+1):
            logger.info(f'Working on IBI iteration {i}/{iterations}')

            with logging_redirect_tqdm(loggers=[logger,devlogger]):
                for p in tqdm(_potentials, desc='Writing CG potentials'):
                    # if 'IBIPotentials/intrabond-1' in p.filename(): import ipdb; ipdb.set_trace()
                    try:    p.write_cg_potential(cg_u, tol=p.tol, box=box)
                    except: p.write_cg_potential(cg_u, box=box)

            if i == 1 and run_minimization:
                logger.info(f'Running brief simulation with small timestep')
                ts0 = engine._get_combined_conf(self).timestep
                engine.simulate( self,
                                 output_name = 'ibi-min', directory = directory,
                                 timestep = ts0/100,
                                 num_steps = 10000, output_period=1000 )

                restart_file = f'{directory}/output/ibi-min.restart'

            name = f'ibi-{i:03d}'
            engine.simulate( self,
                             output_name = name, directory = directory,
                             restart_file = restart_file,
                             replicas = replicas )

            restart_file = f'{directory}/output/{name}{".0" if replicas > 1 else ""}.restart'
            cg_u = _load_cg_u(i)

            with logging_redirect_tqdm(loggers=[logger,devlogger]):
                for p in tqdm(_potentials, desc='Extracting CG distributions'):
                    p.get_cg_distribution(cg_u, box=box, recalculate=False)
                    p.iteration += 1

    def update(self, group , copy=False):
        assert( isinstance(group, Group) )
        if copy:
            group = deepcopy(group)
        group.parent = self
        self.add(group)

    def _get_nonbonded_interaction(self, typeA, typeB):
        for s,A,B in self.nonbonded_interactions:
            if A is None or B is None:
                if A is None and B is None:
                    return s
                elif A is None and typeB.is_same_type(B):
                    return s
                elif B is None and typeA.is_same_type(A):
                    return s
            elif typeA.is_same_type(A) and typeB.is_same_type(B):
                return s

        # raise Exception("No nonbonded scheme found for %s and %s" % (typeA.name, typeB.name))
        # print("WARNING: No nonbonded scheme found for %s and %s" % (typeA.name, typeB.name))

    def _countParticleTypes(self):
        ## TODO: check for modifications to particle that require
        ## automatic generation of new particle type
        type_counts = dict()    # type is key, value is 2-element list of regular particle counts and attached particles

        parts, self.rigid_bodies = [],[]
        for p in self:
            if p.rigid:
                self.rigid_bodies.append(p)
            else:
                parts.append(p)

        for p in parts:
            t = p.type_
            if t in type_counts:
                type_counts[t][0] += 1
            else:
                type_counts[t] = [1,0]

        parts = [p for rb in self if rb.rigid for p in rb.attached_particles]
        for p in parts:
            t = p.type_
            if t in type_counts:
                type_counts[t][1] += 1
            else:
                type_counts[t] = [0,1]

        if len(self.dummy_types) != 0:
            raise("Dummy types have been deprecated")
        # for t in self.dummy_types:
        #     if t not in type_counts:
        #         type_counts[t] = 0

        for i,tA,tB in self.nonbonded_interactions:
            if tA is not None and tA not in type_counts:
                type_counts[tA] = [0,0]
            if tB is not None and tB not in type_counts:
                type_counts[tB] = [0,0]


        self.type_counts = type_counts

        rbtc = dict()
        rbti={}
        for rb in self.rigid_bodies:
            t = rb.type_.name
            if t in rbtc: rbtc[t] += 1
            else:
                rbti[t]=rb.type_  
                rbtc[t] = 1
        self.rigid_body_type_counts = [(k,rbtc[k]) for k in sorted(rbtc.keys())]
        self.rigid_body_index=rbti
        devlogger.debug(f'{self}: Counting types: {type_counts}')
        devlogger.debug(f'{self}: Counting rb types: {rbtc}')

    def _updateParticleOrder(self):
        ## Create ordered list
        self.particles = [p for p in self if not p.rigid]
        self.rigid_bodies = list(sorted([p for p in self if p.rigid], key=lambda x: x.type_))
        # self.particles = sorted(particles, key=lambda p: (p.type_, p.idx))

        ## Update particle indices
        idx = 0
        for p in self.particles:
            p.idx = idx
            idx = idx+1

        ## Add attached particle indices
        # attach particles
        for j,rb in enumerate(self.rigid_bodies):
            for p in rb.attached_particles:
                p.idx = idx
                idx = idx+1

        ## TODO recurse through childrens' group_sites
        for g in self.group_sites:
            g.idx = idx
            idx = idx + 1

        # self.initialCoords = np.array([p.initialPosition for p in self.particles])

    def useNonbondedScheme(self, nbScheme, typeA=None, typeB=None):
        """ deprecated """
        logger.warning('useNonbondedScheme is deprecated! Please update your code to use `add_nonbonded_interaction`')
        self.add_nonbonded_interaction(nbScheme, typeA, typeB)

    def add_nonbonded_interaction(self, nonbonded_potential, typeA=None, typeB=None):
        self.nonbonded_interactions.append( [nonbonded_potential, typeA, typeB] )
        if typeA != typeB:
            self.nonbonded_interactions.append( [nonbonded_potential, typeB, typeA] )

    def prepare_for_simulation(self):
        ...

    def getParticleTypesAndCounts(self):
        """ Includes rigid body-attached particles """
        ## TODO: remove(?)
        if self.type_counts is None:
            self._countParticleTypes()
            self._updateParticleOrder()

        return sorted( self.type_counts.items(), key=lambda x: x[0] )

    def _particleTypePairIter(self):
        typesAndCounts = self.getParticleTypesAndCounts()
        i_skipped = 0
        for i in range(len(typesAndCounts)):
            t1,(n1,rb1) = typesAndCounts[i]
            if n1+rb1 == 0:
                i_skipped += 1
                continue
            j_skipped = 0
            for j in range(i,len(typesAndCounts)):
                t2,(n2,rb2) = typesAndCounts[j]
                if n2+rb2 == 0:
                    j_skipped += 1
                    continue
                if n2 == 0: continue
                yield( [i-i_skipped,j-i_skipped-j_skipped,t1,t2] )

    def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
        raise(NotImplementedError)

    def simulate(self, output_name, **kwargs):
        ## split kwargs
        sim_kws = ['output_directory', 'directory', 'log_file', 'binary', 'num_procs', 'dry_run', 'configuration', 'replicas']
        sim_kwargs = {kw:kwargs[kw] for kw in sim_kws if kw in kwargs}
        engine_kwargs = {k:v for k,v in kwargs.items() if k not in sim_kws}
        engine = ArbdEngine(**engine_kwargs)
        return engine.simulate(self, output_name, **sim_kwargs)

assign_IBI_degrees_of_freedom()

Convenience routine that adds degrees of freedom to corresponding IBI potentials

Source code in mrdna/arbdmodel/submodule/__init__.py
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
def assign_IBI_degrees_of_freedom(self):

    """ Convenience routine that adds degrees of freedom to
    corresponding IBI potentials """

    from .ibi import BondDof, AngleDof, DihedralDof, PairDistributionDof

    self.bonded_ibi_potentials = set()
    self.nonbonded_ibi_potentials = set()

    for get_fn, parts_pot_fn, cls in (
            (self.get_bonds,     lambda x: (x[:2],x[2]), BondDof),
            (self.get_angles,    lambda x: (x[:3],x[3]), AngleDof),
            (self.get_dihedrals, lambda x: (x[:4],x[4]), DihedralDof)
    ):
        for x in get_fn():
            parts, pot = parts_pot_fn(x)
            try: pot.degrees_of_freedom
            except: continue
            pot.degrees_of_freedom.append( cls(*parts) )
            if pot not in self.bonded_ibi_potentials:
                self.bonded_ibi_potentials.add( pot )

    logger.info(f'Gathering nonbonded IBI degrees of freedom')
    all_exclusions = self.get_exclusions()
    try:
        raise               # particle.idx wasn't set without this
        if len(self.type_counts) == 0: raise
    except:
        self.getParticleTypesAndCounts()

    type_to_particles = dict()
    for p in self:
        t = p.type_
        if t not in type_to_particles: type_to_particles[t] = []
        type_to_particles[t].append(p)
    type_to_particles[None] = [p for p in self]

    already_handled = set()
    for pot, tA, tB in self.nonbonded_interactions:
        ## Avoid including particle type pairs that don't apply
        if (tA,tB) in already_handled: continue
        already_handled.add((tA,tB))
        already_handled.add((tB,tA))

        try:
            pot.degrees_of_freedom
        except:
            continue

        def _cond(pair):
            b1 = (tA is None or pair[0].type_ == tA) and (tB is None or pair[1].type_ == tB)
            b2 = (tB is None or pair[0].type_ == tB) and (tA is None or pair[1].type_ == tA)
            return (b1 or b2)

        ex = list(filter(_cond, all_exclusions))
        dof = PairDistributionDof( type_to_particles[tA], type_to_particles[tB], exclusions=ex,
                                   range_=pot.range_, resolution=pot.resolution)

        pot.degrees_of_freedom.append( dof )

        if pot not in self.nonbonded_ibi_potentials:
            self.nonbonded_ibi_potentials.add( pot )

getParticleTypesAndCounts()

Includes rigid body-attached particles

Source code in mrdna/arbdmodel/submodule/__init__.py
1393
1394
1395
1396
1397
1398
1399
1400
def getParticleTypesAndCounts(self):
    """ Includes rigid body-attached particles """
    ## TODO: remove(?)
    if self.type_counts is None:
        self._countParticleTypes()
        self._updateParticleOrder()

    return sorted( self.type_counts.items(), key=lambda x: x[0] )

useNonbondedScheme(nbScheme, typeA=None, typeB=None)

deprecated

Source code in mrdna/arbdmodel/submodule/__init__.py
1380
1381
1382
1383
def useNonbondedScheme(self, nbScheme, typeA=None, typeB=None):
    """ deprecated """
    logger.warning('useNonbondedScheme is deprecated! Please update your code to use `add_nonbonded_interaction`')
    self.add_nonbonded_interaction(nbScheme, typeA, typeB)

ConnectableElement

Abstract base class

Source code in mrdna/arbdmodel/submodule/polymer.py
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
class ConnectableElement():
    """ Abstract base class """
    def __init__(self, connection_locations=None, connections=None):
        if connection_locations is None: connection_locations = []
        if connections is None: connections = []

        ## TODO decide on names
        self.locations = self.connection_locations = connection_locations
        self.connections = connections

    def get_locations(self, type_=None, exclude=()):
        locs = [l for l in self.connection_locations if (type_ is None or l.type_ == type_) and l.type_ not in exclude]
        counter = dict()
        for l in locs:
            if l in counter:
                counter[l] += 1
            else:
                counter[l] = 1
        assert( np.all( [counter[l] == 1 for l in locs] ) )
        return locs

    def get_location_at(self, contour_position, on_fwd_strand=True, new_type="crossover"):
        loc = None
        if (self.num_monomers == 1):
            ## Assumes that intrahelical connections have been made before crossovers
            for l in self.locations:
                if l.on_fwd_strand == on_fwd_strand and l.connection is None:
                    assert(loc is None)
                    loc = l
            # assert( loc is not None )
        else:
            for l in self.locations:
                if l.contour_position == contour_position and l.on_fwd_strand == on_fwd_strand:
                    assert(loc is None)
                    loc = l
        if loc is None:
            loc = Location( self, contour_position=contour_position, type_=new_type, on_fwd_strand=on_fwd_strand )
        return loc

    def get_connections_and_locations(self, connection_type=None, exclude=()):
        """ Returns a list with each entry of the form:
            connection, location_in_self, location_in_other """
        type_ = connection_type
        ret = []
        for c in self.connections:
            if (type_ is None or c.type_ == type_) and c.type_ not in exclude:
                if   c.A.parent is self:
                    ret.append( [c, c.A, c.B] )
                elif c.B.parent is self:
                    ret.append( [c, c.B, c.A] )
                else:
                    import pdb
                    pdb.set_trace()
                    raise Exception("Object contains connection that fails to refer to object")
        return ret

    def _connect(self, other, connection):
        ## TODO fix circular references        
        A,B = [connection.A, connection.B]

        A.connection = B.connection = connection
        self.connections.append(connection)
        if other is not self:
            other.connections.append(connection)
        else:
            raise NotImplementedError("Objects cannot yet be connected to themselves; if you are attempting to make a circular object, try breaking the object into multiple components")
        l = A.parent.locations
        if A not in l: l.append(A)
        l = B.parent.locations
        if B not in l: l.append(B)

get_connections_and_locations(connection_type=None, exclude=())

Returns a list with each entry of the form: connection, location_in_self, location_in_other

Source code in mrdna/arbdmodel/submodule/polymer.py
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
def get_connections_and_locations(self, connection_type=None, exclude=()):
    """ Returns a list with each entry of the form:
        connection, location_in_self, location_in_other """
    type_ = connection_type
    ret = []
    for c in self.connections:
        if (type_ is None or c.type_ == type_) and c.type_ not in exclude:
            if   c.A.parent is self:
                ret.append( [c, c.A, c.B] )
            elif c.B.parent is self:
                ret.append( [c, c.B, c.A] )
            else:
                import pdb
                pdb.set_trace()
                raise Exception("Object contains connection that fails to refer to object")
    return ret

Connection

Class for connecting two elements

Source code in mrdna/arbdmodel/submodule/polymer.py
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
class Connection():
    """ Class for connecting two elements """
    def __init__(self, A, B, type_ = None):
        assert( isinstance(A,Location) )
        assert( isinstance(B,Location) )
        self.A = A
        self.B = B
        self.type_ = type_

    def other(self, location):
        if location is self.A:
            return self.B
        elif location is self.B:
            return self.A
        else:
            raise ValueError

    def delete(self):
        self.A.parent.connections.remove(self)
        if self.B.parent is not self.A.parent:
            self.B.parent.connections.remove(self)
        self.A.connection = None
        self.B.connection = None

    def __repr__(self):
        return "<Connection {}--{}--{}]>".format( self.A, self.type_, self.B )

Location

Site for connection within an object

Source code in mrdna/arbdmodel/submodule/polymer.py
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
class Location():
    """ Site for connection within an object """
    def __init__(self, parent, contour_position, type_, on_fwd_strand = True):
        ## TODO: remove cyclic references(?)
        assert( isinstance(parent, ConnectableElement) )
        self.parent = parent
        self.contour_position = contour_position  # represents position along contour length in polymer
        self.type_ = type_
        # self.particle = None
        self.connection = None

    def get_connected_location(self):
        if self.connection is None:
            return None
        else:
            return self.connection.other(self)

    def set_connection(self, connection):
        self.connection = connection # TODO weakref? 

    def get_monomer_index(self):
        try:
            idx = self.parent.contour_to_monomer_index(self.contour_position, nearest_monomer=True)
        except:
            if self.contour_position == 0:
                idx = 0
            elif self.contour_position == 1:
                idx = self.parent.num_monomers-1
            else:
                raise
        return idx

    def __repr__(self):
        if self.on_fwd_strand:
            on_fwd = "on_fwd_strand"
        else:
            on_fwd = "on_rev_strand"
        return "<Location {}.{}[{:.2f},{:d}]>".format( self.parent.name, self.type_, self.contour_position, self.on_fwd_strand)

PolymerBeads

Bases: Group

Abstract class for bead-based models that are built from Polymer objects

Source code in mrdna/arbdmodel/submodule/polymer.py
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
class PolymerBeads(Group, metaclass=ABCMeta):
    """ Abstract class for bead-based models that are built from Polymer objects """

    def __init__(self, polymer, sequence=None, monomers_per_bead_group = 1, rest_length = None, **kwargs):
        self.polymer = polymer
        self.sequence = sequence
        # self.spring_constant = spring_constant
        self.monomers_per_bead_group = monomers_per_bead_group

        if rest_length is None:
            rest_length = polymer.monomer_length * self.monomers_per_bead_group

        self.rest_length = rest_length

        for prop in ('segname','chain'):
            if prop not in kwargs:
                try:
                    self.__dict__[prop] = polymer.__dict__[prop]
                except:
                    pass

        Group.__init__(self, **kwargs)

    @property
    def monomers_per_bead_group(self):
        return self.__monomers_per_bead_group

    @property
    def num_bead_groups(self):
        return self.__num_bead_groups    

    @monomers_per_bead_group.setter
    def monomers_per_bead_group(self,value):
        if value <= 0:
            raise ValueError("monomers_per_bead_group must be positive")
        self.num_bead_groups = int(np.ceil(self.polymer.num_monomers/value))

    @num_bead_groups.setter
    def num_bead_groups(self,value):
        if value <= 0:
            raise ValueError("num_bead_groups must be positive")
        self.__num_bead_groups = value
        self.__monomers_per_bead_group = self.polymer.num_monomers/value # get exact ratio

    def _clear_beads(self):
        ## self.children = []
        self.clear_all()

    @abstractmethod
    def _generate_ith_bead_group(self, index, position, orientation):
        """ Normally a bead, but sometimes a group of beads """
        bead = PointParticle(type_, position,
                             resid = index+1)
        pass

    @abstractmethod
    def _join_adjacent_bead_groups(self, ids):
        if len(ids) == 2:
            b1,b2 = [self.children[i] for i in ids]
            """ units "10 kJ/N_A" kcal_mol """
            bond = HarmonicBond(k = self.spring_constant,
                                r0 = self.rest_length,
                                rRange = (0,500),
                                resolution = 0.01,
                                maxForce = 50)
            self.add_bond( i=b1, j=b2, bond = bond, exclude=True )
        else:
            pass

    def _generate_beads(self):
        for i in range(self.num_bead_groups):
            c = self.polymer.monomer_index_to_contour(i * self.monomers_per_bead_group + (self.monomers_per_bead_group-1)/2)
            r = self.polymer.contour_to_position(c)
            o = self.polymer.contour_to_orientation(c)

            obj = self._generate_ith_bead_group(i, r, o)
            self.add(obj)

        for di in range(1,4):
            for i in range(len(self.children)-di):
                ids = tuple(i+_di for _di in range(di+1)) 
                self._join_adjacent_bead_groups(ids)

PolymerGroup

Source code in mrdna/arbdmodel/submodule/polymer.py
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
class PolymerGroup():
    def __init__(self, polymers=[],
                 **kwargs):

        """ Simulation system parameters """
        ## TODO decide whether to move these elsewhere
        # self.dimensions = dimensions
        # self.temperature = temperature
        # self.timestep = timestep
        # self.cutoff = cutoff
        # self.decomposition_period = decomposition_period
        # self.pairlist_distance = pairlist_distance
        # self.nonbonded_resolution = nonbonded_resolution
        self.polymers = polymers

        self.grid_potentials = []

        # self.useNonbondedScheme( nbDnaScheme )

    def get_connections(self,type_=None,exclude=()):
        """ Find all connections in model, without double-counting """
        added=set()
        ret=[]
        for s in self.polymers:
            items = [e for e in s.get_connections_and_locations(type_,exclude=exclude) if e[0] not in added]
            added.update([e[0] for e in items])
            ret.extend( list(sorted(items,key=lambda x: x[1].contour_position)) )
        return ret

    def extend(self, other, copy=True, include_strands=False):
        assert( isinstance(other, PolymerSectionModel) )
        if copy:
            for s in other.polymers:
                self.polymers.append(deepcopy(s))
            if include_strands:
                for s in other.strands:
                    self.strands.append(deepcopy(s))
        else:
            for s in other.polymers:
                self.polymers.append(s)
            if include_strands:
                for s in other.strands:
                    self.strands.append(s)
        self._clear_beads()

    def update(self, polymer , copy=False):
        assert( isinstance(polymer, PolymerSection) )
        if copy:
            polymer = deepcopy(polymer)
        self.polymers.append(polymer)

    """ Operations on spline coordinates """
    def translate(self, translation_vector, position_filter=None):
        for s in self.polymers:
            s.translate(translation_vector, position_filter=position_filter)

    def rotate(self, rotation_matrix, about=None, position_filter=None):
        for s in self.polymers:
            s.rotate(rotation_matrix, about=about, position_filter=position_filter)

    def get_center(self, include_ssdna=False):
        if include_ssdna:
            polymers = self.polymers
        else:
            polymers = list(filter(lambda s: isinstance(s,DoubleStrandedPolymerSection), 
                                   self.polymers))
        centers = [s.get_center() for s in polymers]
        weights = [s.num_monomers*2 if isinstance(s,DoubleStrandedPolymerSection) else s.num_monomers for s in polymers]
        # centers,weights = [np.array(a) for a in (centers,weights)]
        return np.average( centers, axis=0, weights=weights)

    def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
        positions = []
        for s in self.polymers:
            positions.append(s.contour_to_position(0))
            positions.append(s.contour_to_position(0.5))
            positions.append(s.contour_to_position(1))
        positions = np.array(positions)
        dx,dy,dz = [(np.max(positions[:,i])-np.min(positions[:,i])+30)*padding_factor for i in range(3)]
        if isotropic:
            dx = dy = dz = max((dx,dy,dz))
        return [dx,dy,dz]

    def add_grid_potential(self, grid_file, scale=1, per_monomer=True):
        grid_file = Path(grid_file)
        if not grid_file.is_file():
            raise ValueError("Grid file {} does not exist".format(grid_file))
        if not grid_file.is_absolute():
            grid_file = Path.cwd() / grid_file
        self.grid_potentials.append((grid_file,scale,per_,pmp,er))

    def vmd_tube_tcl(self, file_name="drawTubes.tcl", radius=5):
        with open(file_name, 'w') as tclFile:
            tclFile.write("## beginning TCL script \n")

            def draw_tube(polymer,radius_value=10, color="cyan", resolution=5):
                tclFile.write("## Tube being drawn... \n")

                contours = np.linspace(0,1, max(2,1+polymer.num_monomers//resolution) )
                rs = [polymer.contour_to_position(c) for c in contours]

                radius_value = str(radius_value)
                tclFile.write("graphics top color {} \n".format(str(color)))
                for i in range(len(rs)-2):
                    r0 = rs[i]
                    r1 = rs[i+1]
                    filled = "yes" if i in (0,len(rs)-2) else "no"
                    tclFile.write("graphics top cylinder {{ {} {} {} }} {{ {} {} {} }} radius {} resolution 30 filled {} \n".format(r0[0], r0[1], r0[2], r1[0], r1[1], r1[2], str(radius_value), filled))
                    tclFile.write("graphics top sphere {{ {} {} {} }} radius {} resolution 30\n".format(r1[0], r1[1], r1[2], str(radius_value)))
                r0 = rs[-2]
                r0 = rs[-1]
                tclFile.write("graphics top cylinder {{ {} {} {} }} {{ {} {} {} }} radius {} resolution 30 filled yes \n".format(r0[0], r0[1], r0[2], r1[0], r1[1], r1[2], str(radius_value)))

            ## material
            tclFile.write("graphics top materials on \n")
            tclFile.write("graphics top material AOEdgy \n")

            ## iterate through the model polymers
            for s in self.polymers:
                draw_tube(s,radius,"cyan")


    def vmd_cylinder_tcl(self, file_name="drawCylinders.tcl", radius=5):

        with open(file_name, 'w') as tclFile:
            tclFile.write("## beginning TCL script \n")
            def draw_cylinder(polymer,radius_value=10,color="cyan"):
                tclFile.write("## cylinder being drawn... \n")
                r0 = polymer.contour_to_position(0)
                r1 = polymer.contour_to_position(1)

                radius_value = str(radius_value)
                color = str(color)

                tclFile.write("graphics top color {} \n".format(color))
                tclFile.write("graphics top cylinder {{ {} {} {} }} {{ {} {} {} }} radius {} resolution 30 filled yes \n".format(r0[0], r0[1], r0[2], r1[0], r1[1], r1[2], radius_value))

            ## material
            tclFile.write("graphics top materials on \n")
            tclFile.write("graphics top material AOEdgy \n")

            ## iterate through the model polymers
            for s in self.polymers:
                draw_cylinder(s,radius,"cyan")

__init__(polymers=[], **kwargs)

Simulation system parameters

Source code in mrdna/arbdmodel/submodule/polymer.py
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
def __init__(self, polymers=[],
             **kwargs):

    """ Simulation system parameters """
    ## TODO decide whether to move these elsewhere
    # self.dimensions = dimensions
    # self.temperature = temperature
    # self.timestep = timestep
    # self.cutoff = cutoff
    # self.decomposition_period = decomposition_period
    # self.pairlist_distance = pairlist_distance
    # self.nonbonded_resolution = nonbonded_resolution
    self.polymers = polymers

    self.grid_potentials = []

get_connections(type_=None, exclude=())

Find all connections in model, without double-counting

Source code in mrdna/arbdmodel/submodule/polymer.py
425
426
427
428
429
430
431
432
433
def get_connections(self,type_=None,exclude=()):
    """ Find all connections in model, without double-counting """
    added=set()
    ret=[]
    for s in self.polymers:
        items = [e for e in s.get_connections_and_locations(type_,exclude=exclude) if e[0] not in added]
        added.update([e[0] for e in items])
        ret.extend( list(sorted(items,key=lambda x: x[1].contour_position)) )
    return ret

PolymerModel

Bases: ArbdModel

Source code in mrdna/arbdmodel/submodule/polymer.py
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
class PolymerModel(ArbdModel, metaclass=ABCMeta):
    def __init__(self, polymers,
                 sequences = None,
                 monomers_per_bead_group = 1,
                 **kwargs):

        """ Assign sequences """
        if sequences is None:
            sequences = [None for i in range(len(polymers))]

        self.polymer_group = PolymerGroup(polymers)
        self.sequences = sequences
        self.monomers_per_bead_group = monomers_per_bead_group
        ArbdModel.__init__(self, [], **kwargs)

        """ Generate beads """
        self.generate_beads()

    def update_splines(self, coords):
        i = 0
        for p,c in zip(self.polymer_group.polymers,self.children):
            n = len(c.children)
            p.set_splines(np.linspace(0,1,n), coords[i:i+n])
            i += n

    @abstractmethod
    def _generate_polymer_beads(self, polymer, sequence, polymer_index = None):
        pass

    def generate_beads(self):
        self.clear_all()
        self.peptides = [self._generate_polymer_beads(p,s,i)
                         for i,(p,s) in enumerate(zip(self.polymer_group.polymers,self.sequences))]

        for s in self.peptides:
            self.add(s)
            s._generate_beads()

__init__(polymers, sequences=None, monomers_per_bead_group=1, **kwargs)

Assign sequences

Source code in mrdna/arbdmodel/submodule/polymer.py
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
def __init__(self, polymers,
             sequences = None,
             monomers_per_bead_group = 1,
             **kwargs):

    """ Assign sequences """
    if sequences is None:
        sequences = [None for i in range(len(polymers))]

    self.polymer_group = PolymerGroup(polymers)
    self.sequences = sequences
    self.monomers_per_bead_group = monomers_per_bead_group
    ArbdModel.__init__(self, [], **kwargs)

    """ Generate beads """
    self.generate_beads()

PolymerSection

Bases: ConnectableElement

Base class that describes a linear section of a polymer

Source code in mrdna/arbdmodel/submodule/polymer.py
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
class PolymerSection(ConnectableElement):
    """ Base class that describes a linear section of a polymer """

    def __init__(self, name, num_monomers,
                 monomer_length,
                 start_position = None,
                 end_position = None, 
                 parent = None,
                 **kwargs):

        ConnectableElement.__init__(self, connection_locations=[], connections=[])

        if 'segname' not in kwargs:
            self.segname = name

        for key,val in kwargs.items():
            self.__dict__[key] = val

        self.num_monomers = int(num_monomers)
        self.monomer_length = monomer_length

        if start_position is None: start_position = np.array((0,0,0))

        if end_position is None:
            end_position = np.array((0,0,self.monomer_length*num_monomers)) + start_position
        self.start_position = start_position
        self.end_position = end_position

        self.start_orientation = None

        ## Set up interpolation for positions
        self._set_splines_from_ends()

        # self.sequence = None

    def __repr__(self):
        return "<{} {}[{:d}]>".format( type(self), self.segname, self.num_monomers )

    def set_splines(self, contours, coords):
        tck, u = interpolate.splprep( coords.T, u=contours, s=0, k=1)
        self.position_spline_params = (tck,u)

    def set_orientation_splines(self, contours, quaternions):
        tck, u = interpolate.splprep( quaternions.T, u=contours, s=0, k=1)
        self.quaternion_spline_params = (tck,u)

    def get_center(self):
        tck, u = self.position_spline_params
        return np.mean(self.contour_to_position(u), axis=0)

    def _get_location_positions(self):
        return [self.contour_to_monomer_index(l.contour_position) for l in self.locations]

    def insert_monomers(self, at_monomer: int, num_monomer: int, seq=tuple()):
        assert(np.isclose(np.around(num_monomer),num_monomer))
        if at_monomer < 0:
            raise ValueError("Attempted to insert DNA into {} at a negative location".format(self))
        if at_monomer > self.num_monomer-1:
            raise ValueError("Attempted to insert DNA into {} at beyond the end of the Polymer".format(self))
        if num_monomers < 0:
            raise ValueError("Attempted to insert DNA a negative amount of DNA into {}".format(self))

        num_monomers = np.around(num_monomer)
        monomer_positions = self._get_location_positions()
        new_monomer_positions = [p if p <= at_monomer else p+num_monomers for p in monomer_positions]

        ## TODO: handle sequence

        self.num_monomers = self.num_monomer+num_monomer

        for l,p in zip(self.locations, new_monomer_positions):
            l.contour_position = self.monomer_index_to_contour(p)

    def remove_monomers(self, first_monomer: int, last_monomer: int):
        """ Removes nucleotides between first_monomer and last_monomer, inclusive """
        assert(np.isclose(np.around(first_monomer),first_monomer))
        assert(np.isclose(np.around(last_monomer),last_monomer))
        tmp = min((first_monomer,last_monomer))
        last_monomer = max((first_monomer,last_monomer))
        fist_monomer = tmp

        if first_monomer < 0 or first_monomer > self.num_monomer-2:
            raise ValueError("Attempted to remove DNA from {} starting at an invalid location {}".format(self, first_monomer))
        if last_monomer < 1 or last_monomer > self.num_monomer-1:
            raise ValueError("Attempted to remove DNA from {} ending at an invalid location {}".format(self, last_monomer))
        if first_monomer == last_monomer:
            return

        first_monomer = np.around(first_monomer)
        last_monomer = np.around(last_monomer)

        monomer_positions = self._get_location_positions()

        bad_locations = list(filter(lambda p: p >= first_monomer and p <= last_monomer, monomer_positions))
        if len(bad_locations) > 0:
            raise Exception("Attempted to remove DNA containing locations {} from {} between {} and {}".format(bad_locations,self,first_monomer,last_monomer))

        removed_monomer = last_monomer-first_monomer+1
        new_monomer_positions = [p if p <= last_monomer else p-removed_monomer for p in monomer_positions]
        num_monomers = self.num_monomer-removed_monomer

        if self.sequence is not None and len(self.sequence) == self.num_monomer:
            self.sequence = [s for s,i in zip(self.sequence,range(self.num_monomer)) 
                                if i < first_monomer or i > last_monomer]
            assert( len(self.sequence) == num_monomers )

        self.num_monomers = num_monomers

        for l,p in zip(self.locations, new_monomer_positions):
            l.contour_position = self.monomer_position_to_contour(p)

    def __filter_contours(contours, positions, position_filter, contour_filter):
        u = contours
        r = positions

        ## Filter
        ids = list(range(len(u)))
        if contour_filter is not None:
            ids = list(filter(lambda i: contour_filter(u[i]), ids))
        if position_filter is not None:
            ids = list(filter(lambda i: position_filter(r[i,:]), ids))
        return ids

    def translate(self, translation_vector, position_filter=None, contour_filter=None):
        dr = np.array(translation_vector)
        tck, u = self.position_spline_params
        r = self.contour_to_position(u)

        ids = PolymerSection.__filter_contours(u, r, position_filter, contour_filter)
        if len(ids) == 0: return

        ## Translate
        r[ids,:] = r[ids,:] + dr[np.newaxis,:]
        self.set_splines(u,r)

    def rotate(self, rotation_matrix, about=None, position_filter=None, contour_filter=None):
        tck, u = self.position_spline_params
        r = self.contour_to_position(u)

        ids = PolymerSection.__filter_contours(u, r, position_filter, contour_filter)
        if len(ids) == 0: return

        if about is None:
            ## TODO: do this more efficiently
            r[ids,:] = np.array([rotation_matrix.dot(r[i,:]) for i in ids])
        else:
            dr = np.array(about)
            ## TODO: do this more efficiently
            r[ids,:] = np.array([rotation_matrix.dot(r[i,:]-dr) + dr for i in ids])

        self.set_splines(u,r)

        if self.quaternion_spline_params is not None:
            ## TODO: performance: don't shift between quaternion and matrix representations so much
            tck, u = self.quaternion_spline_params
            orientations = [self.contour_to_orientation(v) for v in u]
            for i in ids:
                orientations[i,:] = rotation_matrix.dot(orientations[i])
            quats = [quaternion_from_matrix(o) for o in orientations]
            self.set_orientation_splines(u, quats)

    def _set_splines_from_ends(self, resolution=4):
        self.quaternion_spline_params = None
        r0 = np.array(self.start_position)[np.newaxis,:]
        r1 = np.array(self.end_position)[np.newaxis,:]
        u = np.linspace(0,1, max(3,self.num_monomers//int(resolution)))
        s = u[:,np.newaxis]
        coords = (1-s)*r0 + s*r1
        self.set_splines(u, coords)

    def contour_to_monomer_index(self, contour_pos, nearest_monomer=False):
        idx = contour_pos*(self.num_monomer) - 0.5
        if nearest_monomer:
            assert( np.isclose(np.around(idx),idx) )
            idx = np.around(idx)
        return idx

    def monomer_index_to_contour(self,monomer_pos):
        return (monomer_pos+0.5)/(self.num_monomers)

    def contour_to_position(self,s):
        p = interpolate.splev( s, self.position_spline_params[0] )
        if len(p) > 1: p = np.array(p).T
        return p

    def contour_to_tangent(self,s):
        t = interpolate.splev( s, self.position_spline_params[0], der=1 )
        t = (t / np.linalg.norm(t,axis=0))
        return t.T


    def contour_to_orientation(self,s):
        assert( isinstance(s,float) or isinstance(s,int) or len(s) == 1 )   # TODO make vectorized version

        if self.quaternion_spline_params is None:
            axis = self.contour_to_tangent(s)
            axis = axis / np.linalg.norm(axis)
            rotAxis = np.cross(axis,np.array((0,0,1)))
            rotAxisL = np.linalg.norm(rotAxis)
            zAxis = np.array((0,0,1))

            if rotAxisL > 0.001:
                theta = np.arcsin(rotAxisL) * 180/np.pi
                if axis.dot(zAxis) < 0: theta = 180-theta
                orientation0 = rotationAboutAxis( rotAxis/rotAxisL, theta, normalizeAxis=False ).T
            else:
                orientation0 = np.eye(3) if axis.dot(zAxis) > 0 else \
                               rotationAboutAxis( np.array((1,0,0)), 180, normalizeAxis=False )
            if self.start_orientation is not None:
                orientation0 = orientation0.dot(self.start_orientation)

            # orientation = rotationAboutAxis( axis, self.twist_per_nt*self.contour_to_monomer_index(s), normalizeAxis=False )
            # orientation = orientation.dot(orientation0)
            orientation = orientation0
        else:
            q = interpolate.splev( s, self.quaternion_spline_params[0] )
            if len(q) > 1: q = np.array(q).T # TODO: is this needed?
            orientation = quaternion_to_matrix(q)

        return orientation

    def get_contour_sorted_connections_and_locations(self,type_):
        sort_fn = lambda c: c[1].contour_position
        cl = self.get_connections_and_locations(type_)
        return sorted(cl, key=sort_fn)

    def add_location(self, nt, type_, on_fwd_strand=True):
        ## Create location if needed, add to polymer
        c = self.monomer_index_to_contour(nt)
        assert(c >= 0 and c <= 1)
        # TODO? loc = self.Location( contour_position=c, type_=type_, on_fwd_strand=is_fwd )
        loc = Location( self, contour_position=c, type_=type_, on_fwd_strand=on_fwd_strand )
        self.locations.append(loc)

    def iterate_connections_and_locations(self, reverse=False):
        ## connections to other polymers
        cl = self.get_contour_sorted_connections_and_locations()
        if reverse:
            cl = cl[::-1]

        for c in cl:
            yield c

remove_monomers(first_monomer, last_monomer)

Removes nucleotides between first_monomer and last_monomer, inclusive

Source code in mrdna/arbdmodel/submodule/polymer.py
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def remove_monomers(self, first_monomer: int, last_monomer: int):
    """ Removes nucleotides between first_monomer and last_monomer, inclusive """
    assert(np.isclose(np.around(first_monomer),first_monomer))
    assert(np.isclose(np.around(last_monomer),last_monomer))
    tmp = min((first_monomer,last_monomer))
    last_monomer = max((first_monomer,last_monomer))
    fist_monomer = tmp

    if first_monomer < 0 or first_monomer > self.num_monomer-2:
        raise ValueError("Attempted to remove DNA from {} starting at an invalid location {}".format(self, first_monomer))
    if last_monomer < 1 or last_monomer > self.num_monomer-1:
        raise ValueError("Attempted to remove DNA from {} ending at an invalid location {}".format(self, last_monomer))
    if first_monomer == last_monomer:
        return

    first_monomer = np.around(first_monomer)
    last_monomer = np.around(last_monomer)

    monomer_positions = self._get_location_positions()

    bad_locations = list(filter(lambda p: p >= first_monomer and p <= last_monomer, monomer_positions))
    if len(bad_locations) > 0:
        raise Exception("Attempted to remove DNA containing locations {} from {} between {} and {}".format(bad_locations,self,first_monomer,last_monomer))

    removed_monomer = last_monomer-first_monomer+1
    new_monomer_positions = [p if p <= last_monomer else p-removed_monomer for p in monomer_positions]
    num_monomers = self.num_monomer-removed_monomer

    if self.sequence is not None and len(self.sequence) == self.num_monomer:
        self.sequence = [s for s,i in zip(self.sequence,range(self.num_monomer)) 
                            if i < first_monomer or i > last_monomer]
        assert( len(self.sequence) == num_monomers )

    self.num_monomers = num_monomers

    for l,p in zip(self.locations, new_monomer_positions):
        l.contour_position = self.monomer_position_to_contour(p)