segmentmodel.py 67 KB
Newer Older
1001
1002
1003
1004
        self.children = self.segments = segments

        self._bonded_potential = dict() # cache bonded potentials

1005
        self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
1006
1007
1008

        self.useNonbondedScheme( nbDnaScheme )

1009

cmaffeo2's avatar
cmaffeo2 committed
1010
    def get_connections(self,type_=None,exclude=[]):
1011
1012
1013
1014
        """ Find all connections in model, without double-counting """
        added=set()
        ret=[]
        for s in self.segments:
cmaffeo2's avatar
cmaffeo2 committed
1015
            items = [e for e in s.get_connections_and_locations(type_,exclude=exclude) if e[0] not in added]
1016
1017
1018
            added.update([e[0] for e in items])
            ret.extend( items )
        return ret
1019
    
1020
    def _recursively_get_beads_within_bonds(self,b1,bonds,done=[]):
1021
        ret = []
1022
1023
1024
1025
1026
1027
1028
1029
1030
        done = list(done)
        done.append(b1)
        if bonds == 0:
            return [[]]

        for b2 in b1.intrahelical_neighbors:
            if b2 in done: continue
            for tmp in self._recursively_get_beads_within_bonds(b2, bonds-1, done):
                ret.append( [b2]+tmp )
1031
1032
        return ret

1033
1034
1035
1036
1037
    def _get_intrahelical_beads(self,num=2):
        ## TODO: add check that this is not called before adding intrahelical_neighbors in _generate_bead_model

        assert(num >= 2)

1038
1039
        ret = []
        for s in self.segments:
1040
1041
1042
1043
1044
            for b1 in s.beads:
                for bead_list in self._recursively_get_beads_within_bonds(b1, num-1):
                    assert(len(bead_list) == num-1)
                    if b1.idx < bead_list[-1].idx: # avoid double-counting
                        ret.append([b1]+bead_list)
1045
1046
        return ret

1047

1048
1049
    def _get_intrahelical_angle_beads(self):
        return self._get_intrahelical_beads(num=3)
1050

cmaffeo2's avatar
cmaffeo2 committed
1051
    def _get_potential(self, type_, kSpring, d, max_potential = None):
1052
        key = (type_, kSpring, d, max_potential)
1053
1054
        if key not in self._bonded_potential:
            if type_ == "bond":
1055
                self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,500), max_potential=max_potential)
1056
            elif type_ == "angle":
cmaffeo2's avatar
cmaffeo2 committed
1057
1058
                self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
                # , resolution = 1, maxForce=0.1)
1059
            elif type_ == "dihedral":
cmaffeo2's avatar
cmaffeo2 committed
1060
                self._bonded_potential[key] = HarmonicDihedral(kSpring,d, max_potential=max_potential)
1061
1062
1063
1064
1065
1066
1067
            else:
                raise Exception("Unhandled potential type '%s'" % type_)
        return self._bonded_potential[key]
    def get_bond_potential(self, kSpring, d):
        return self._get_potential("bond", kSpring, d)
    def get_angle_potential(self, kSpring, d):
        return self._get_potential("angle", kSpring, d)
cmaffeo2's avatar
cmaffeo2 committed
1068
1069
1070
1071
    def get_dihedral_potential(self, kSpring, d, max_potential=None):
        while d > 180: d-=360
        while d < -180: d+=360
        return self._get_potential("dihedral", kSpring, d, max_potential)
1072
1073


cmaffeo2's avatar
cmaffeo2 committed
1074
1075
    def _getParent(self, *beads ):
        if np.all( [b1.parent == b2.parent 
1076
                    for b1,b2 in zip(beads[:-1],beads[1:])] ):
cmaffeo2's avatar
cmaffeo2 committed
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
            return beads[0].parent
        else:
            return self

    def _get_twist_spring_constant(self, sep):
        """ sep in nt """
        kT = 0.58622522         # kcal/mol
        twist_persistence_length = 90  # set semi-arbitrarily as there is a large spread in literature
        ## <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
        ##   Assume A is small
        ## int[B_] :=  Normal[Integrate[ Series[Cos[x] Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]/
        ##             Integrate[Series[Exp[-B x^2], {B, 0, 1}], {x, 0, \[Pi]}]]

        ## Actually, without assumptions I get fitFun below
        ## From http://www.annualreviews.org/doi/pdf/10.1146/annurev.bb.17.060188.001405
        ##   units "3e-19 erg cm/ 295 k K" "nm" =~ 73
        Lp = twist_persistence_length/0.34 

        fitFun = lambda x: np.real(erf( (4*np.pi*x + 1j)/(2*np.sqrt(x)) )) * np.exp(-1/(4*x)) / erf(2*np.sqrt(x)*np.pi) - np.exp(-sep/Lp)
        k = opt.leastsq( fitFun, x0=np.exp(-sep/Lp) )
        return k[0][0] * 2*kT*0.00030461742
1098
1099
1100
1101
1102
1103

    """ Mapping between different resolution models """
    def _clear_beads(self):
        for s in self.segments:
            s.clear_all()
        self.clear_all(keep_children=True)
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
        assert( len([b for b in self]) == 0 )
        locParticles = []
        # for c,A,B in self.get_connections():
        for s in self.segments:
            for c,A,B in s.get_connections_and_locations():
                for l in (A,B):
                    if l.particle is not None:
                        locParticles.append(A.particle)
        assert( len(locParticles) == 0 )
        assert( len([b for s in self.segments for b in s.beads]) == 0 )
1114
1115
1116
1117

    def _update_segment_positions(self, bead_coordinates):
        """ Set new function for each segments functions
        contour_to_position and contour_to_orientation """
cmaffeo2's avatar
cmaffeo2 committed
1118

1119
        for s in self.segments:
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
            cabs = s.get_connections_and_locations("intrahelical")
            beads = list(set(s.beads + [A.particle for c,A,B in cabs]))

            ## Add nearby beads
            for c,A,B in cabs:
                ## TODOTODO test?
                bs = filter( B.particle.intrahelical_neighbors, lambda x: x is not None and x not in beads )
                beads.extend(bs)
                bs = filter( bs.intrahelical_neighbors, lambda x: x is not None and x not in beads )
                beads.extend(bs)

                        
1132
1133
            ## Skip beads that are None (some locations were not assigned a particle to avoid double-counting) 
            beads = [b for b in beads if b is not None]
1134
1135
1136
1137
1138
1139
            contours = [b.get_contour_position(s) for b in beads]

            cb = sorted( zip(contours,beads), key=lambda a:a[0] )
            beads = [b for c,b in cb] 
            contours = [c for c,b in cb] 

1140
1141
1142
1143
            ids = [b.idx for b in beads]
            
            """ Get positions """
            positions = bead_coordinates[ids,:].T
1144
1145
1146
1147
1148
1149

            try:
                tck, u = interpolate.splprep( positions, u=contours, s=0, k=3 )
            except:
                tck, u = interpolate.splprep( positions, u=contours, s=0, k=1 )

1150
1151
1152
            s.position_spline_params = tck

            """ Get twist """
1153
1154
1155
1156
1157
1158
            cb = [e for e in cb if 'orientation_bead' in e[1].__dict__]
            beads = [b for c,b in cb] 
            contours = [c for c,b in cb] 
            ids = [b.idx for b in beads]
            # if 'orientation_bead' in beads[0].__dict__:
            if len(beads) > 3:
1159
1160
                tangents = s.contour_to_tangent(contours)
                quats = []
cmaffeo2's avatar
cmaffeo2 committed
1161
                lastq = None
1162
1163
                for b,t in zip(beads,tangents):
                    o = b.orientation_bead
1164
1165
1166
                    positions
                    # angleVec = o.position - b.position
                    angleVec = bead_coordinates[o.idx,:] - bead_coordinates[b.idx,:]
1167
1168
1169
                    angleVec = angleVec - angleVec.dot(t)*t
                    angleVec = angleVec/np.linalg.norm(angleVec)
                    y = np.cross(t,angleVec)
cmaffeo2's avatar
cmaffeo2 committed
1170
1171
                    assert( np.abs(np.linalg.norm(y) - 1) < 1e-2 )
                    q = quaternion_from_matrix( np.array([angleVec,y,t]).T)
1172
1173
                    # q = quaternion_from_matrix( np.array([angleVec,y,t]))
                    
cmaffeo2's avatar
cmaffeo2 committed
1174
1175
1176
1177
                    if lastq is not None:
                        if q.dot(lastq) < 0:
                            q = -q
                    quats.append( q )
1178
                    lastq = q
cmaffeo2's avatar
cmaffeo2 committed
1179

1180
                # pdb.set_trace()
1181
                quats = np.array(quats)
1182
                ## TODOTODO test smoothing
1183
                try:
1184
                    tck, u = interpolate.splprep( quats.T, u=contours, s=3, k=3 )
1185
1186
                except:
                    tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=1 )
1187
1188
1189
                s.quaternion_spline_params = tck

    def _generate_bead_model(self,
1190
                             max_basepairs_per_bead = 7,
cmaffeo2's avatar
cmaffeo2 committed
1191
                             max_nucleotides_per_bead = 4,
1192
1193
1194
                             local_twist=False,
                             escapable_twist=True):

1195

1196
        segments = self.segments
1197
1198
        for s in segments:
            s.local_twist = local_twist
1199
1200

        """ Simplify connections """
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
        # d_nt = dict()           # 
        # for s in segments:
        #     d_nt[s] = 1.5/(s.num_nts-1)
        # for s in segments:
        #     ## replace consecutive crossovers with
        #     cl = sorted( s.get_connections_and_locations("crossover"), key=lambda x: x[1].address )
        #     last = None
        #     for entry in cl:
        #         c,A,B = entry
        #         if last is not None and \
        #            (A.address - last[1].address) < d_nt[s]:
        #             same_type = c.type_ == last[0].type_
        #             same_dest_seg = B.container == last[2].container
        #             if same_type and same_dest_seg:
        #                 if np.abs(B.address - last[2].address) < d_nt[B.container]:
        #                     ## combine
        #                     A.combine = last[1]
        #                     B.combine = last[2]

        #                     ...
        #         # if last is not None:
        #         #     s.bead_locations.append(last)
        #         ...
        #         last = entry
        # del d_nt
1226

1227
        """ Generate beads at intrahelical junctions """
1228
        if self.DEBUG: print( "Adding intrahelical beads at junctions" )
1229

1230
        ## Loop through all connections, generating beads at appropriate locations
1231
1232
1233
        for c,A,B in self.get_connections("intrahelical"):
            s1,s2 = [l.container for l in (A,B)]

1234
1235
1236
1237
            ## TODO be more elegant!
            if A.on_fwd_strand == False: continue # TODO verify this avoids double-counting
            assert( A.particle is None )
            assert( B.particle is None )
1238

1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
            ## TODO: offload the work here to s1
            a1,a2 = [l.address   for l in (A,B)]            
            # a1,a2 = [a - s.nt_pos_to_contour(0.5) if a == 0 else a + s.nt_pos_to_contour(0.5) for a,s in zip((a1,a2),(s1,s2))]
            for a in (a1,a2):
                assert( a in (0,1,0.0,1.0) )
            a1,a2 = [a - s.nt_pos_to_contour(0.5) if a == 0 else a + s.nt_pos_to_contour(0.5) for a,s in zip((a1,a2),(s1,s2))]
            
            ## TODO improve this for combinations of ssDNA and dsDNA (maybe a1/a2 should be calculated differently)
            if isinstance(s2,DoubleStrandedSegment):
                b = s2._generate_one_bead(a2,0)
1249
            else:
1250
1251
                b = s1._generate_one_bead(a1,0)
            A.particle = B.particle = b
1252

1253
        """ Generate beads at other junctions """
cmaffeo2's avatar
cmaffeo2 committed
1254
1255
1256
1257
1258
1259
1260
        for c,A,B in self.get_connections(exclude="intrahelical"):
            s1,s2 = [l.container for l in (A,B)]
            if A.particle is not None and B.particle is not None:
                continue
            assert( A.particle is None )
            assert( B.particle is None )

1261
            ## TODO: offload the work here to s1/s2 (?)
cmaffeo2's avatar
cmaffeo2 committed
1262
            a1,a2 = [l.address   for l in (A,B)]
1263
1264

            b = s1.get_nearest_bead(a1)
1265
            if b is not None and s1.contour_to_nt_pos(np.abs(b.contour_position-a1)) < 1.9:
1266
1267
1268
1269
1270
1271
1272
                ## combine beads
                b.contour_position = 0.5*(b.contour_position + a1) # avg position
                A.particle = b
            else:
                A.particle = s1._generate_one_bead(a1,0)

            b = s2.get_nearest_bead(a2)
1273
            if b is not None and s2.contour_to_nt_pos(np.abs(b.contour_position-a2)) < 1.9:
1274
1275
1276
1277
1278
                ## combine beads
                b.contour_position = 0.5*(b.contour_position + a2) # avg position
                B.particle = b
            else:
                B.particle = s2._generate_one_bead(a2,0)
cmaffeo2's avatar
cmaffeo2 committed
1279

1280
1281
1282
1283
1284
1285
1286
        """ Some tests """
        for c,A,B in self.get_connections("intrahelical"):
            for l in (A,B):
                if l.particle is None: continue
                assert( l.particle.parent is not None )

        """ Generate beads in between """
1287
        if self.DEBUG: print("Generating beads")
1288
        for s in segments:
cmaffeo2's avatar
cmaffeo2 committed
1289
            s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
1290

1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
        # """ Combine beads at junctions as needed """
        # for c,A,B in self.get_connections():
        #    ...

        """ Add intrahelical neighbors at connections """
        for c,A,B in self.get_connections("intrahelical"):
            b1,b2 = [l.particle for l in (A,B)]
            if b1 is b2:
                ## already handled by Segment._generate_beads
                continue
            else:
                for b in (b1,b2): assert( b is not None )
1303
1304
1305
1306
1307
                if b2 not in b1.intrahelical_neighbors:
                    b1.intrahelical_neighbors.append(b2)
                    b2.intrahelical_neighbors.append(b1)
                assert(len(b1.intrahelical_neighbors) <= 2)
                assert(len(b2.intrahelical_neighbors) <= 2)
1308
1309

        """ Reassign bead types """
1310
        if self.DEBUG: print("Assigning bead types")
1311
1312
1313
        beadtype_s = dict()
        for segment in segments:
            for b in segment:
1314
                b.num_nts = np.around( b.num_nts, decimals=1 )
1315
1316
1317
1318
1319
1320
1321
1322
1323
                key = (b.type_.name[0].upper(), b.num_nts)
                if key in beadtype_s:
                    b.type_ = beadtype_s[key]
                else:
                    t = deepcopy(b.type_)
                    if key[0] == "D":
                        t.__dict__["nts"] = b.num_nts*2
                    elif key[0] == "S":
                        t.__dict__["nts"] = b.num_nts
cmaffeo2's avatar
cmaffeo2 committed
1324
1325
                    elif key[0] == "O":
                        t.__dict__["nts"] = b.num_nts
1326
1327
                    else:
                        raise Exception("TODO")
cmaffeo2's avatar
cmaffeo2 committed
1328
                    # print(t.nts)
1329
                    t.name = t.name + "%03d" % (10*t.nts)
1330
1331
                    beadtype_s[key] = b.type_ = t

1332
1333
1334
        """ Update bead indices """
        self._countParticleTypes() # probably not needed here
        self._updateParticleOrder()
1335

1336
        """ Add intrahelical bond potentials """
1337
        if self.DEBUG: print("Adding intrahelical bond potentials")
cmaffeo2's avatar
cmaffeo2 committed
1338
        dists = dict()          # intrahelical distances built for later use
1339
1340
1341
        intra_beads = self._get_intrahelical_beads() 
        if self.DEBUG: print("  Adding %d bonds" % len(intra_beads))
        for b1,b2 in intra_beads:
cmaffeo2's avatar
cmaffeo2 committed
1342
            parent = self._getParent(b1,b2)
1343
1344
1345

            ## TODO: could be sligtly smarter about sep
            sep = 0.5*(b1.num_nts+b2.num_nts)
cmaffeo2's avatar
cmaffeo2 committed
1346

1347
            conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
1348
            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D":
1349
                elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
1350
                d = 3.4*sep
1351
                k = conversion*elastic_modulus/d
1352
            else:
1353
1354
                ## TODO: get better numbers our ssDNA model
                elastic_modulus = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
1355
                d = 5*sep
1356
                k = conversion*elastic_modulus/d
1357
                # print(sep,d,k)
1358
1359
              
            if b1 not in dists:
cmaffeo2's avatar
cmaffeo2 committed
1360
                dists[b1] = dict()
1361
            if b2 not in dists:
cmaffeo2's avatar
cmaffeo2 committed
1362
1363
1364
1365
1366
                dists[b2] = dict()
            # dists[b1].append([b2,sep])
            # dists[b2].append([b1,sep])
            dists[b1][b2] = sep
            dists[b2][b1] = sep
1367
1368

            # dists[[b1,b2]] = dists[[b2,b1]] = sep
1369
1370
            bond = self.get_bond_potential(k,d)
            parent.add_bond( b1, b2, bond, exclude=True )
1371
1372

        """ Add intrahelical angle potentials """
1373
        if self.DEBUG: print("Adding intrahelical angle potentials")
1374
        for b1,b2,b3 in self._get_intrahelical_angle_beads():
1375
1376
            ## TODO: could be slightly smarter about sep
            sep = 0.5*b1.num_nts+b2.num_nts+0.5*b3.num_nts
cmaffeo2's avatar
cmaffeo2 committed
1377
            parent = self._getParent(b1,b2,b3)
1378

1379
1380
1381
1382
1383
1384
            kT = 0.58622522         # kcal/mol
            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D" and b3.type_.name[0] == "D":
                ## <cos(q)> = exp(-s/Lp) = integrate( x^4 exp(-A x^2) / 2, {x, 0, pi} ) / integrate( x^2 exp(-A x^2), {x, 0, pi} )
                ## <cos(q)> ~ 1 - 3/4A                                                                                            
                ## where A = k_spring / (2 kT)                                                                                    
                k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
cmaffeo2's avatar
cmaffeo2 committed
1385
                if local_twist:
1386
                    ## TODO optimize this paramter
cmaffeo2's avatar
cmaffeo2 committed
1387
                    k *= 0.5    # halve because orientation beads have similar springs
1388
1389
            else:
                ## TODO: get correct number from ssDNA model
cmaffeo2's avatar
cmaffeo2 committed
1390
                k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2
1391

1392
1393
1394
1395
            angle = self.get_angle_potential(k,180)
            parent.add_angle( b1, b2, b3, angle )

        """ Add intrahelical exclusions """
1396
        if self.DEBUG: print("Adding intrahelical exclusions")
1397
1398
1399
        beads = dists.keys()
        def _recursively_get_beads_within(b1,d,done=[]):
            ret = []
cmaffeo2's avatar
cmaffeo2 committed
1400
            for b2,sep in dists[b1].items():
1401
1402
1403
1404
1405
1406
1407
                if b2 in done: continue
                if sep < d:
                    ret.append( b2 )
                    done.append( b2 )
                    tmp = _recursively_get_beads_within(b2, d-sep, done)
                    if len(tmp) > 0: ret.extend(tmp)
            return ret
1408

1409
1410
1411
        exclusions = set()
        for b1 in beads:
            exclusions.update( [(b1,b) for b in _recursively_get_beads_within(b1, 20, done=[b1])] )
1412
1413

        if self.DEBUG: print("Adding %d exclusions" % len(exclusions))
1414
        for b1,b2 in exclusions:
cmaffeo2's avatar
cmaffeo2 committed
1415
            parent = self._getParent(b1,b2)
1416
            parent.add_exclusion( b1, b2 )
cmaffeo2's avatar
cmaffeo2 committed
1417
1418
1419

        """ Twist potentials """
        if local_twist:
1420
            if self.DEBUG: print("Adding twist potentials")
cmaffeo2's avatar
cmaffeo2 committed
1421
1422
1423

            for b1 in beads:
                if "orientation_bead" not in b1.__dict__: continue
cmaffeo2's avatar
cmaffeo2 committed
1424
                for b2,sep in dists[b1].items():
cmaffeo2's avatar
cmaffeo2 committed
1425
                    if "orientation_bead" not in b2.__dict__: continue
1426
                    if b2.idx < b1.idx: continue # Don't double-count
cmaffeo2's avatar
cmaffeo2 committed
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445

                    p1,p2 = [b.parent for b in (b1,b2)]
                    o1,o2 = [b.orientation_bead for b in (b1,b2)]

                    parent = self._getParent( b1, b2 )

                    """ Add heuristic 90 degree potential to keep orientation bead orthogonal """
                    k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
                    pot = self.get_angle_potential(k,90)
                    parent.add_angle(o1,b1,b2, pot)
                    parent.add_angle(b1,b2,o2, pot)
                                        
                    ## TODO: improve this
                    twist_per_nt = 0.5 * (p1.twist_per_nt + p2.twist_per_nt)
                    angle = sep*twist_per_nt
                    if angle > 360 or angle < -360:
                        raise Exception("The twist between beads is too large")
                        
                    k = self._get_twist_spring_constant(sep)
1446
1447
1448
1449
                    if escapable_twist:
                        pot = self.get_dihedral_potential(k,angle,max_potential=1)
                    else:
                        pot = self.get_dihedral_potential(k,angle)
cmaffeo2's avatar
cmaffeo2 committed
1450
1451
                    parent.add_dihedral(o1,b1,b2,o2, pot)

1452
        """ Add connection potentials """
1453
        for c,A,B in self.get_connections("terminal_crossover"):
cmaffeo2's avatar
cmaffeo2 committed
1454
            ## TODO: use a better description here
1455
1456
1457
            b1,b2 = [loc.particle for loc in (c.A,c.B)]
            pot = self.get_bond_potential(4,18.5)
            self.add_bond(b1,b2, pot)
1458

1459
        
1460
        crossover_bead_pots = set()
cmaffeo2's avatar
cmaffeo2 committed
1461
1462
        for c,A,B in self.get_connections("crossover"):
            b1,b2 = [loc.particle for loc in (c.A,c.B)]
1463
1464
1465
1466
1467
1468

            ## Avoid double-counting
            if (b1,b2) in crossover_bead_pots: continue
            crossover_bead_pots.add((b1,b2))
            crossover_bead_pots.add((b2,b1))
                
cmaffeo2's avatar
cmaffeo2 committed
1469
1470
1471
1472
            pot = self.get_bond_potential(4,18.5)
            self.add_bond(b1,b2, pot)

            if not local_twist:
1473
1474
                ## TODOTODO
                
cmaffeo2's avatar
cmaffeo2 committed
1475
1476
1477
1478
                ...
            else:
                k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
                if 'orientation_bead' in b1.__dict__:
1479
1480
1481
                    # t0 = 90 + 60
                    t0 = 150
                    if A.on_fwd_strand: t0 = 30 # TODO handle antiparallel segments
cmaffeo2's avatar
cmaffeo2 committed
1482
1483
1484
1485
                    o = b1.orientation_bead
                    pot = self.get_angle_potential(k,t0)
                    self.add_angle( o,b1,b2, pot )
                else:
1486
1487
                    t0 = 150
                    if B.on_fwd_strand: t0 = 30
cmaffeo2's avatar
cmaffeo2 committed
1488
1489
1490
1491
                    o = b2.orientation_bead
                    pot = self.get_angle_potential(k,t0)
                    self.add_angle( b1,b2,o, pot )

1492
1493
1494
                ## Get beads above and below
                u1,u2 = [b.get_intrahelical_above() for b in (b1,b2)]
                d1,d2 = [b.get_intrahelical_below() for b in (b1,b2)]
cmaffeo2's avatar
cmaffeo2 committed
1495

cmaffeo2's avatar
cmaffeo2 committed
1496
                k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
cmaffeo2's avatar
cmaffeo2 committed
1497
1498
1499
1500
                t0 = 90
                if 'orientation_bead' in b1.__dict__:
                    o1 = b1.orientation_bead
                    if u2 is not None:
cmaffeo2's avatar
cmaffeo2 committed
1501
1502
                        k = k_fn( dists[b2][u2] )
                        pot = self.get_dihedral_potential(k,t0)
cmaffeo2's avatar
cmaffeo2 committed
1503
1504
                        self.add_dihedral( o1,b1,b2,u2, pot )
                    elif d2 is not None:
cmaffeo2's avatar
cmaffeo2 committed
1505
1506
1507
                        k = k_fn( dists[b2][d2] )
                        pot = self.get_dihedral_potential(k,t0)
                        self.add_dihedral( o1,b1,b2,d2, pot )
cmaffeo2's avatar
cmaffeo2 committed
1508
1509
1510
                if 'orientation_bead' in b2.__dict__:
                    o2 = b2.orientation_bead
                    if u1 is not None:
cmaffeo2's avatar
cmaffeo2 committed
1511
1512
                        k = k_fn( dists[b1][u1] )
                        pot = self.get_dihedral_potential(k,t0)
cmaffeo2's avatar
cmaffeo2 committed
1513
1514
                        self.add_dihedral( o2,b2,b1,u1, pot )
                    elif d1 is not None:
cmaffeo2's avatar
cmaffeo2 committed
1515
1516
                        k = k_fn( dists[b1][d1] )
                        pot = self.get_dihedral_potential(k,t0)
cmaffeo2's avatar
cmaffeo2 committed
1517
1518
1519
1520
1521
                        self.add_dihedral( o2,b2,b1,d1, pot )

                if 'orientation_bead' in b1.__dict__ and 'orientation_bead' in b2.__dict__:
                    if u1 is not None and u2 is not None:
                        t0 = 0
cmaffeo2's avatar
cmaffeo2 committed
1522
                        k = k_fn( 0.5*(dists[b1][u1]+dists[b2][u2]) )
cmaffeo2's avatar
cmaffeo2 committed
1523
1524
1525
1526
                        pot = self.get_dihedral_potential(k,t0)
                        self.add_dihedral( u1,b1,b2,u2,  pot )
                    elif d1 is not None and d2 is not None:
                        t0 = 0
cmaffeo2's avatar
cmaffeo2 committed
1527
                        k = k_fn( 0.5*(dists[b1][d1]+dists[b2][d2]) )
cmaffeo2's avatar
cmaffeo2 committed
1528
1529
1530
1531
                        pot = self.get_dihedral_potential(k,t0)
                        self.add_dihedral( d1,b1,b2,d2, pot )
                    elif d1 is not None and u2 is not None:
                        t0 = 180
cmaffeo2's avatar
cmaffeo2 committed
1532
                        k = k_fn( 0.5*(dists[b1][d1]+dists[b2][u2]) )
cmaffeo2's avatar
cmaffeo2 committed
1533
1534
1535
1536
                        pot = self.get_dihedral_potential(k,t0)
                        self.add_dihedral( d1,b1,b2,u2, pot )
                    elif u1 is not None and d2 is not None:
                        t0 = 180
cmaffeo2's avatar
cmaffeo2 committed
1537
                        k = k_fn( 0.5*(dists[b1][u1]+dists[b2][d2]) )
cmaffeo2's avatar
cmaffeo2 committed
1538
1539
                        pot = self.get_dihedral_potential(k,t0)
                        self.add_dihedral( u1,b1,b2,d2, pot )
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551

        def get_consecutive_crossovers():
            processed_segments = set()
            for s in self.segments:
                if s in processed_segments: continue
            
        for s in self.segments:
            for cl in s.get_contour_sorted_connections_and_locations("crossover"):
                
            
        

1552
    def _generate_strands(self):
1553
        self.strands = strands = []
cmaffeo2's avatar
cmaffeo2 committed
1554

1555
1556
1557
1558
1559
        """ Ensure unconnected ends have 5prime Location objects """
        for seg in self.segments:
            ## TODO move into Segment calls
            import pdb
            if seg.start5.connection is None:
1560
1561
1562
1563
1564
1565
1566
                add_end = True
                for l in seg.get_locations("5prime"):
                    if l.address == 0 and l.on_fwd_strand:
                        add_end = False
                        break
                if add_end:
                    seg.add_5prime(0) 
1567
            if 'end5' in seg.__dict__ and seg.end5.connection is None:
1568
1569
1570
1571
1572
1573
1574
                add_end = True
                for l in seg.get_locations("5prime"):
                    if l.address == 1 and (l.on_fwd_strand is False):
                        add_end = False
                        break
                if add_end:
                    seg.add_5prime(seg.num_nts-1,on_fwd_strand=False)
1575
1576

            if 'start3' in seg.__dict__ and seg.start3.connection is None:
1577
1578
1579
1580
1581
1582
1583
                add_end = True
                for l in seg.get_locations("3prime"):
                    if l.address == 0 and (l.on_fwd_strand is False):
                        add_end = False
                        break
                if add_end:
                    seg.add_3prime(0,on_fwd_strand=False)
1584
            if seg.end3.connection is None:
1585
1586
1587
1588
1589
1590
1591
                add_end = True
                for l in seg.get_locations("3prime"):
                    if l.address == 1 and l.on_fwd_strand:
                        add_end = False
                        break
                if add_end:
                    seg.add_3prime(seg.num_nts-1)
1592
1593
1594
1595
1596
1597
1598
1599

            # print( [(l,l.get_connected_location()) for l in seg.locations] )
            # addresses = np.array([l.address for l in seg.get_locations("5prime")])
            # if not np.any( addresses == 0 ):
            #     ## check if end is connected
        # for c,l,B in self.get_connections_and_locations():
        #     if c[0]

cmaffeo2's avatar
cmaffeo2 committed
1600
        """ Build strands from connectivity of helices """
1601
1602
1603
1604
1605
        def _recursively_build_strand(strand, segment, pos, is_fwd, mycounter=0):
            mycounter+=1
            if mycounter > 1000:
                import pdb
                pdb.set_trace()
cmaffeo2's avatar
cmaffeo2 committed
1606
            s,seg = [strand, segment]
1607
1608

            end_pos, next_seg, next_pos, next_dir = seg.get_strand_segment(pos, is_fwd)
cmaffeo2's avatar
cmaffeo2 committed
1609
            s.add_dna(seg, pos, end_pos, is_fwd)
1610

cmaffeo2's avatar
cmaffeo2 committed
1611
1612
            if next_seg is not None:
                # print("  next_dir: {}".format(next_dir))
1613
                _recursively_build_strand(s, next_seg, next_pos, next_dir, mycounter)
cmaffeo2's avatar
cmaffeo2 committed
1614
1615
1616

        for seg in self.segments:
            locs = seg.get_5prime_locations()
1617
1618
1619
            if locs is None: continue
            # for pos, is_fwd in locs:
            for l in locs:
1620
1621
                # print("Tracing",l)
                pos = seg.contour_to_nt_pos(l.address)
1622
                is_fwd = l.on_fwd_strand
cmaffeo2's avatar
cmaffeo2 committed
1623
1624
                s = Strand()
                _recursively_build_strand(s, seg, pos, is_fwd)
1625
                # print("{} {}".format(seg.name,s.num_nts))
cmaffeo2's avatar
cmaffeo2 committed
1626
1627
                strands.append(s)
        
1628
        self.strands = sorted(strands, key=lambda s:s.num_nts)[::-1] # or something
1629
1630
1631
1632
1633
1634
1635
        ## relabel segname
        counter = 0
        for s in self.strands:
            if s.segname is None:
                s.segname = "D%03d" % counter
            counter += 1

1636
1637
1638
1639
    def _update_orientations(self,orientation):
        for s in self.strands:
            s.update_atomic_orientations(orientation)

1640

1641
    def _generate_atomic_model(self):
1642
        self.children = self.strands
1643
        for s in self.strands:
cmaffeo2's avatar
cmaffeo2 committed
1644
            s.generate_atomic_model()
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
        return

        ## Angle optimization
        angles = np.linspace(-180,180,180)
        score = []
        for a in angles:
            o = rotationAboutAxis([0,0,1], a)
            sum2 = count = 0
            for s in self.strands:
                s.update_atomic_orientations(o)
                for s1,s2 in zip(s.strand_segments[:-1],s.strand_segments[1:]):
                    nt1 = s1.children[-1]
                    nt2 = s2.children[0]
                    o3 = nt1.atoms_by_name["O3'"]
                    p = nt2.atoms_by_name["P"]
                    sum2 += np.sum((p.collapsedPosition()-o3.collapsedPosition())**2)
                    count += 1
            score.append(sum2/count)
        print(angles[np.argmin(score)])
        print(score)
For faster browsing, not all history is shown. View entire blame