segmentmodel.py 25.8 KB
Newer Older
1
2
import numpy as np
from arbdmodel import PointParticle, ParticleType, Group, ArbdModel
cmaffeo2's avatar
cmaffeo2 committed
3
from coords import rotationAboutAxis
4
5
6
7
from nonbonded import *
from copy import copy, deepcopy
from nbPot import nbDnaScheme

cmaffeo2's avatar
cmaffeo2 committed
8
9
10
from scipy.special import erf
import scipy.optimize as opt

11
12

"""
cmaffeo2's avatar
cmaffeo2 committed
13
14
15
TODO:
 - document
 - 
16
17
18
19
20
21
22
23
24
25
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


"""

class Location():
    """ Site for connection within an object """
    def __init__(self, container, address, type_):
        self.container = container
        self.address = address
        self.type_ = type_
        self.particle = None

class Connection():
    """ Abstract base class for connection between 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_
        
# class ConnectableElement(Transformable):
class ConnectableElement():
    """ Abstract base class """
    def __init__(self, connections=[]):
        self.connections = connections

    def _connect(self, other, connection):
        self.connections.append(connection)
        other.connections.append(connection)

    def _find_connections(self, loc):
        return [c for c in self.connections if c.A == loc or c.B == loc]
       
class Segment(ConnectableElement, Group):

    """ Base class that describes a segment of DNA. When built from
    cadnano models, should not span helices """

    """Define basic particle types"""
    dsDNA_particle = ParticleType("D",
                                  diffusivity = 43.5,
                                  mass = 300,
                                  radius = 3,                 
                              )
cmaffeo2's avatar
cmaffeo2 committed
61
62
63
64
65
66
67
    orientation_particle = ParticleType("O",
                                        diffusivity = 100,
                                        mass = 300,
                                        radius = 1,
                                    )
    # orientation_bond = HarmonicBond(10,2)
    orientation_bond = HarmonicBond(30,1.5)
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

    ssDNA_particle = ParticleType("S",
                                  diffusivity = 43.5,
                                  mass = 150,
                                  radius = 3,                 
                              )

    def __init__(self, name, num_nts, 
                 start_position = np.array((0,0,0)),
                 end_position = None, 
                 segment_model = None):

        Group.__init__(self, name, children=[])
        ConnectableElement.__init__(self, connections=[])

cmaffeo2's avatar
cmaffeo2 committed
83
84
85
86
87
        self.start_orientation = None
        self.twist_per_nt = 0

        self.beads = [c for c in self.children] # self.beads will not contain orientation beads

88
89
90
91
92
93
94
95
96
        self._bead_model_generation = 0    # TODO: remove?
        self.segment_model = segment_model # TODO: remove?

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

cmaffeo2's avatar
cmaffeo2 committed
97
98
99
100
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead ):
        raise NotImplementedError

    def _generate_one_bead(self, pos, nts, orientation=None):
101
102
103
104
105
106
107
        raise NotImplementedError

    def _assign_particles_to_locations(self):
        raise NotImplementedError

    def get_all_consecutive_beads(self, number):
        assert(number >= 1)
cmaffeo2's avatar
cmaffeo2 committed
108
        ## Assume that consecutive beads in self.beads are bonded
109
        ret = []
cmaffeo2's avatar
cmaffeo2 committed
110
111
        for i in range(len(self.beads)-number+1):
            tmp = [self.beads[i+j] for j in range(0,number)]
112
113
114
            ret.append( tmp )
        return ret

115
    def get_beads_before_bead(self, bead, number, inclusive=False):
cmaffeo2's avatar
cmaffeo2 committed
116
117
118
        ## Assume that consecutive beads in self.beads are bonded
        i = self.beads.index(bead)
        l = len(self.beads)
119
120
121
122
123
        if i-number < 0:
            raise Exception("Not enough beads after bead")
        
        start = 1
        if inclusive: start = 0
cmaffeo2's avatar
cmaffeo2 committed
124
        return [self.beads[i-j] for j in range(start,number)]
125

126
    def get_beads_after_bead(self, bead, number, inclusive=False):
cmaffeo2's avatar
cmaffeo2 committed
127
128
129
        ## Assume that consecutive beads in self.beads are bonded
        i = self.beads.index(bead)
        l = len(self.beads)
130
131
132
133
134
        if i+number >= l:
            raise Exception("Not enough beads after bead")
        
        start = 1
        if inclusive: start = 0
cmaffeo2's avatar
cmaffeo2 committed
135
        return [self.beads[i+i] for j in range(start,number)]
136

137
138
139
140
    # def get_bead_pairs_within(self, cutoff):
    #     for b1,b2 in self.get_all_consecutive_beads(self, number)


cmaffeo2's avatar
cmaffeo2 committed
141
    def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
142
143
144
        
        """ Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """

cmaffeo2's avatar
cmaffeo2 committed
145
146
        ## TODO: decide whether to remove bead_model argument

147
148
149
150
        # self._bead_model_generation += 1
        # self._bead_model_max_nts_per_bead = max_nts_per_bead

        direction = self.end_position - self.start_position
cmaffeo2's avatar
cmaffeo2 committed
151
        num_beads = self._get_num_beads( max_basepairs_per_bead, max_nucleotides_per_bead )
152
        nts_per_bead = float(self.num_nts)/num_beads
cmaffeo2's avatar
cmaffeo2 committed
153
        twist_per_bead = nts_per_bead * self.twist_per_nt
154
155
156
157
158
159
160
161
162
163

        last = None

        for i in range(num_beads+1):
            nts = nts_per_bead
            if i == 0 or i == num_beads: 
                nts *= 0.5

            s = i*float(nts_per_bead)/(self.num_nts) # contour
            pos = direction * s + self.start_position
cmaffeo2's avatar
cmaffeo2 committed
164
165
166
167
168
169
170
171
            if self.start_orientation is not None:
                axis = self.start_orientation.dot( np.array((0,0,1)) )
                orientation = rotationAboutAxis( axis, s*self.twist_per_nt*self.num_nts, normalizeAxis=False )
                ## TODO: ensure this is correct
                # orientation = self.start_orientation.dot(orientation) # .dot( self.start_orientation )
                orientation = orientation.dot( self.start_orientation )
            else:
                orientation = None
172

cmaffeo2's avatar
cmaffeo2 committed
173
            b = self._generate_one_bead(pos,nts,orientation)
174
            self.children.append(b)
cmaffeo2's avatar
cmaffeo2 committed
175
176
177
178
179
180
181
            self.beads.append(b) # don't add orientation bead

            if "orientation_bead" in b.__dict__:
                o = b.orientation_bead
                self.children.append(o)
                self.add_bond(b,o, Segment.orientation_bond)
                
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
            # if last is not None:
            #     self.add_bond( i=last, j=b, bond="ssdna" )
            # last = b
        self._assign_particles_to_locations()

    def _regenerate_beads(self, max_nts_per_bead=4, ):
        ...

    def _generate_atomic(self, atomic_model):
        ...
    

class DoubleStrandedSegment(Segment):

    """ Class that describes a segment of ssDNA. When built from
    cadnano models, should not span helices """

    def __init__(self, name, num_nts, start_position = np.array((0,0,0)),
                 end_position = None, 
                 segment_model = None,
cmaffeo2's avatar
cmaffeo2 committed
202
203
                 local_twist = False,
                 num_turns = None,
204
                 start_orientation = None):
cmaffeo2's avatar
cmaffeo2 committed
205
206
207
        
        self.helical_rise = 10.44
        self.distance_per_nt = 3.4
208
209
210
211
212
        Segment.__init__(self, name, num_nts, 
                         start_position,
                         end_position, 
                         segment_model)

cmaffeo2's avatar
cmaffeo2 committed
213
214
215
216
217
218
219
220
221
        self.local_twist = local_twist
        if num_turns is None:
            num_turns = float(num_nts) / self.helical_rise
        self.twist_per_nt = float(360 * num_turns) / num_nts

        if start_orientation is None:
            start_orientation = np.array(((1,0,0),(0,1,0),(0,0,1)))
        self.start_orientation = start_orientation

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
        self.nicks = []

        self.start5 = Location( self, address=0, type_= "end5" )
        self.start3 = Location( self, address=0, type_ = "end3" )

        self.end5 = Location( self, address=-1, type_= "end5" )
        self.end3 = Location( self, address=-1, type_ = "end3" )

    ## Convenience methods
    def connect_start5(self, end3, force_connection=False):
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
        self._connect_ends( self.start5, end3, force_connection = force_connection )
    def connect_start3(self, end5, force_connection=False):
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
        self._connect_ends( self.start3, end5, force_connection = force_connection )
    def connect_end3(self, end5, force_connection=False):
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
        self._connect_ends( self.end3, end5, force_connection = force_connection )
    def connect_end5(self, end3, force_connection=False):
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
        self._connect_ends( self.end5, end3, force_connection = force_connection )

        
    ## Real work
    def _connect_ends(self, end1, end2, force_connection):

        ## validate the input
        for end in (end1, end2):
            assert( isinstance(end, Location) )
            assert( end.type_ in ("end3","end5") )
        assert( end1.type_ != end2.type_ )

        end1.container._connect( end2.container, Connection( end1, end2, type_="intrahelical" ) )

cmaffeo2's avatar
cmaffeo2 committed
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
        return (self.num_nts // max_basepairs_per_bead) + 1

    def _generate_one_bead(self, pos, nts, orientation=None):
        if self.local_twist:
            assert(orientation is not None)
            # opos = pos + np.array((2,0,0)).dot(orientation)
            opos = pos + orientation.dot( np.array((Segment.orientation_bond.r0,0,0)) )
            o = PointParticle( Segment.orientation_particle, opos, nts,
                               num_nts=nts, parent=self )
            bead = PointParticle( Segment.dsDNA_particle, pos, nts,
                                  num_nts=nts, parent=self, orientation_bead=o )

        else:
            bead = PointParticle( Segment.dsDNA_particle, pos, nts,
                                  num_nts=nts, parent=self )
        return bead

278
279

    def _assign_particles_to_locations(self):
cmaffeo2's avatar
cmaffeo2 committed
280
281
        self.start3.particle =  self.start5.particle = self.beads[0]
        self.end3.particle   =  self.end5.particle   = self.beads[-1]
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

    def _generate_atomic(self, atomic_model):
        ...
    
        
    # def add_crossover(self, locationInA, B, locationInB):
    #     j = Crossover( [self, B], [locationInA, locationInB] )
    #     self._join(B,j)

    # def add_internal_crossover(self, locationInA, B, locationInB):
    #     j = Crossover( [self, B], [locationInA, locationInB] )
    #     self._join(B,j)


    # def stack_end(self, myEnd):
    #     ## Perhaps this should not really be possible; these ends should be part of same helix
    #     ...

    # def connect_strand(self, other):
    #     ...
        
    # def break_apart(self):
    #     """Break into smaller pieces so that "crossovers" are only at the ends"""
    #     ...

class SingleStrandedSegment(Segment):

    """ Class that describes a segment of ssDNA. When built from
    cadnano models, should not span helices """

    def __init__(self, name, num_nts, start_position = np.array((0,0,0)),
                 end_position = None, 
                 segment_model = None):

        self.distance_per_nt = 5
        Segment.__init__(self, name, num_nts, 
                         start_position,
                         end_position, 
                         segment_model)

        self.start = self.end5 = Location( self, address=0, type_= "end5" )
        self.end = self.end3 = Location( self, address=-1, type_ = "end3" )

    def connect_3end(self, end5, force_connection=False):
        self._connect_end( end5,  _5_to_3 = False, force_connection = force_connection )

    def connect_5end(self, end3, force_connection=False):
        self._connect_end( end3,  _5_to_3 = True, force_connection = force_connection )

    def _connect_end(self, other, _5_to_3, force_connection):
        assert( isinstance(other, Location) )
        if _5_to_3 == True:
            my_end = self.end5
            assert( other.type_ == "end3" )
        else:
            my_end = self.end3
            assert( other.type_ == "end5" )

        self._connect( other.container, Connection( my_end, other, type_="intrahelical" ) )

cmaffeo2's avatar
cmaffeo2 committed
342
343
344
345
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
        return (self.num_nts // max_nucleotides_per_bead) + 1

    def _generate_one_bead(self, pos, nts, orientation=None):
346
347
348
349
350
351
352
353
354
355
356
357
        return PointParticle( Segment.ssDNA_particle, pos, nts,
                              num_nts=nts, parent=self )

    def _assign_particles_to_locations(self):
        self.start.particle = self.children[0]
        self.end.particle = self.children[-1]

    def _generate_atomic(self, atomic_model):
        ...
    

class SegmentModel(ArbdModel):
cmaffeo2's avatar
cmaffeo2 committed
358
359
360
    def __init__(self, segments=[], local_twist=True,
                 max_basepairs_per_bead=7,
                 max_nucleotides_per_bead=4,
361
362
363
364
365
366
367
368
369
370
371
                 dimensions=(1000,1000,1000), temperature=291,
                 timestep=50e-6, cutoff=50, 
                 decompPeriod=10000, pairlistDistance=None, 
                 nonbondedResolution=0):


        ArbdModel.__init__(self,segments,
                           dimensions, temperature, timestep, cutoff, 
                           decompPeriod, pairlistDistance=None,
                           nonbondedResolution=0)

cmaffeo2's avatar
cmaffeo2 committed
372
373
        # self.max_basepairs_per_bead = max_basepairs_per_bead     # dsDNA
        # self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
374
375
376
377
        self.children = self.segments = segments

        self._bonded_potential = dict() # cache bonded potentials

cmaffeo2's avatar
cmaffeo2 committed
378
        self._generate_bead_model(segments, max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
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
405
    
    def _get_intrahelical_beads(self):
        ret = []
        for s in self.segments:
            ret.extend( s.get_all_consecutive_beads(2) )

        for s in self.segments:
            for c in s.connections:
                if c.type_ == "intrahelical":
                    if c.A.container == s: # avoid double-counting
                        b1,b2 = [loc.particle for loc in (c.A,c.B)]
                        for b in (b1,b2): assert( b is not None )
                        ret.append( [b1,b2] )
        return ret

    def _get_intrahelical_angle_beads(self):
        ret = []
        for s in self.segments:
            ret.extend( s.get_all_consecutive_beads(3) )

        for s1 in self.segments:
            for c in s1.connections:
                if c.A.container != s1: continue
                s2 = c.B.container
                if c.type_ == "intrahelical":
                    b1,b2 = [loc.particle for loc in (c.A,c.B)]
                    for b in (b1,b2): assert( b is not None )
406
407
                    ## TODO: make this code more robust
                    ## TODO: find bug that makes things be out of order
408
                    try:
409
410
411
                        b0 = s1.get_beads_before_bead(b1,1)
                        assert(len(b0) == 1)
                        b0 = b0[0]
412
413
414
415
416
                        assert( b0 is not None )
                        ret.append( [b0,b1,b2] )
                    except:
                        ...
                    try:
417
418
419
                        b0 = s1.get_beads_after_bead(b1,1)
                        assert(len(b0) == 1)
                        b0 = b0[0]
420
421
422
423
424
                        assert( b0 is not None )
                        ret.append( [b2,b1,b0] )
                    except:
                        ...
                    try:
425
426
427
                        b3 = s2.get_beads_before_bead(b2,1)
                        assert(len(b3) == 1)
                        b3 = b3[0]
428
429
430
431
432
                        assert( b3 is not None )
                        ret.append( [b3,b2,b1] )
                    except:
                        ...
                    try:
433
434
435
                        b3 = s2.get_beads_after_bead(b2,1)
                        assert(len(b3) == 1)
                        b3 = b3[0]
436
437
438
439
440
441
                        assert( b3 is not None )
                        ret.append( [b1,b2,b3] )
                    except:
                        ...
        return ret

442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
    # def _get_intrahelical_bead_pairs_within(self, cutoff):
        
    #     dist = dict()
    #     for b1,b2 in self._get_intrahelical_beads:
    #         dist(b1,b2)


    #     ret = []
    #     for s in self.segments:
    #         ret.extend( s.get_bead_pairs_within(cutoff) )

    #     for s1 in self.segments:
    #         for c in s1.connections:
    #             if c.A.container != s1: continue
    #             s2 = c.B.container
    #             if c.type_ == "intrahelical":

    #     ret

cmaffeo2's avatar
cmaffeo2 committed
461
    def _get_potential(self, type_, kSpring, d, max_potential = None):
462
463
464
        key = (type_,kSpring,d)
        if key not in self._bonded_potential:
            if type_ == "bond":
cmaffeo2's avatar
cmaffeo2 committed
465
                self._bonded_potential[key] = HarmonicBond(kSpring,d, max_potential=max_potential)
466
            elif type_ == "angle":
cmaffeo2's avatar
cmaffeo2 committed
467
468
                self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
                # , resolution = 1, maxForce=0.1)
469
            elif type_ == "dihedral":
cmaffeo2's avatar
cmaffeo2 committed
470
                self._bonded_potential[key] = HarmonicDihedral(kSpring,d, max_potential=max_potential)
471
472
473
474
475
476
477
            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
478
479
480
481
    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)
482
483


cmaffeo2's avatar
cmaffeo2 committed
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
    def _getParent(self, *beads ):
        ## TODO: test
        if np.all( [b1.parent == b2.parent 
                    for b1,b2 in zip(beads[::2],beads[1::2])] ):
            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
        

511
512
    def _generate_bead_model(self, segments,
                             max_basepairs_per_bead = 7,
cmaffeo2's avatar
cmaffeo2 committed
513
514
515
                             max_nucleotides_per_bead = 4,
                             local_twist=False
    ):
516
517
518
519


        """ Generate beads """
        for s in segments:
cmaffeo2's avatar
cmaffeo2 committed
520
521
522
            if local_twist:
                s.local_twist = True
            s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542

        """ Combine beads at junctions as needed """
        for s in segments:
            for c in s.connections:
                if c.A.container == s:
                    ...

        """ Reassign bead types """
        beadtype_s = dict()
        for segment in segments:
            for b in segment:
                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
543
544
                    elif key[0] == "O":
                        t.__dict__["nts"] = b.num_nts
545
546
                    else:
                        raise Exception("TODO")
cmaffeo2's avatar
cmaffeo2 committed
547
                    # print(t.nts)
548
549
550
551
                    t.name = t.name + "%03d" % (100*t.nts)
                    beadtype_s[key] = b.type_ = t


552
553
        """ Add intrahelical bond potentials """
        dists = dict()          # built for later use
554
        for b1,b2 in self._get_intrahelical_beads():
cmaffeo2's avatar
cmaffeo2 committed
555
            parent = self._getParent(b1,b2)
556
557
558
559
            if b1.parent == b2.parent:
                sep = 0.5*(b1.num_nts+b2.num_nts)
            else:
                sep = 1
cmaffeo2's avatar
cmaffeo2 committed
560

561
            conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
562
            if b1.type_.name[0] == "D" and b1.type_.name[0] == "D":
563
                elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
564
                d = 3.4*sep
565
                k = conversion*elastic_modulus/d
566
            else:
567
568
                ## TODO: get better numbers our ssDNA model
                elastic_modulus = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
569
                d = 5*sep
570
571
572
573
574
575
576
577
578
579
                k = conversion*elastic_modulus/d
              
            if b1 not in dists:
                dists[b1] = []
            if b2 not in dists:
                dists[b2] = []
            dists[b1].append([b2,sep])
            dists[b2].append([b1,sep])

            # dists[[b1,b2]] = dists[[b2,b1]] = sep
580
581
            bond = self.get_bond_potential(k,d)
            parent.add_bond( b1, b2, bond, exclude=True )
582
583
584

        """ Add intrahelical angle potentials """
        for b1,b2,b3 in self._get_intrahelical_angle_beads():
585
            
586
587
588
589
590
591
592
593
594
595
            sep = 0
            if b1.parent == b2.parent:
                sep += 0.5*(b1.num_nts+b2.num_nts)
            else:
                sep += 1
            if b2.parent == b3.parent:
                sep += 0.5*(b2.num_nts+b3.num_nts)
            else:
                sep += 1
                
cmaffeo2's avatar
cmaffeo2 committed
596
            parent = self._getParent(b1,b2,b3)
597

598
599
600
601
602
603
            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
604
605
                if local_twist:
                    k *= 0.5    # halve because orientation beads have similar springs
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
            else:
                ## TODO: get correct number from ssDNA model
                k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2                
            # k *= 1e-6
            angle = self.get_angle_potential(k,180)
            parent.add_angle( b1, b2, b3, angle )

        """ Add intrahelical exclusions """
        beads = dists.keys()
        def _recursively_get_beads_within(b1,d,done=[]):
            ret = []
            for b2,sep in dists[b1]:
                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
            
        exclusions = set()
        for b1 in beads:
            exclusions.update( [(b1,b) for b in _recursively_get_beads_within(b1, 20, done=[b1])] )
        for b1,b2 in exclusions:
cmaffeo2's avatar
cmaffeo2 committed
630
            parent = self._getParent(b1,b2)
631
            parent.add_exclusion( b1, b2 )
cmaffeo2's avatar
cmaffeo2 committed
632
633
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
671

        """ Twist potentials """
        if local_twist:

            ## TODO: decide whether to add bond here
            # """Add bonds between orientation bead and parent"""
            # for s in self.segments:
            #     for b,o in zip(s.children[::2],s.children[1::2]):
            #         s.add_bond(
                    

            for b1 in beads:
                if "orientation_bead" not in b1.__dict__: continue
                for b2,sep in dists[b1]:
                    if "orientation_bead" not in b2.__dict__: continue

                    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)
                    pot = self.get_dihedral_potential(k,angle,max_potential=1)
                    parent.add_dihedral(o1,b1,b2,o2, pot)

                
                
672
        
673
674
675
676
677
678
679
680
681
682
683
684
        """ Add connection potentials """
        ## TODO


    # def get_bead(self, location):
    #     if type(location.container) is not list:
    #         s = self.segments.index(location.container)
    #         s.get_bead(location.address)
    #     else:
    #         r
    #         ...