segmentmodel.py 27.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
    orientation_particle = ParticleType("O",
                                        diffusivity = 100,
                                        mass = 300,
                                        radius = 1,
                                    )
    # orientation_bond = HarmonicBond(10,2)
67
    orientation_bond = HarmonicBond(30,1.5, rRange = (0,500) )
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
        self._bead_model_generation = 0    # TODO: remove?
        self.segment_model = segment_model # TODO: remove?

cmaffeo2's avatar
cmaffeo2 committed
91
        self.num_nts = int(num_nts)
92
93
94
95
96
        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
151
152
        

cmaffeo2's avatar
cmaffeo2 committed
153
        num_beads = self._get_num_beads( max_basepairs_per_bead, max_nucleotides_per_bead )
154
        nts_per_bead = float(self.num_nts)/num_beads
cmaffeo2's avatar
cmaffeo2 committed
155
        twist_per_bead = nts_per_bead * self.twist_per_nt
156
157
158

        last = None

159
160
161
162
163
164
165
166
167
168
169
170
171
        if num_beads <= 2:
            ## not yet implemented for dsDNA
            assert( isinstance(self, SingleStrandedSegment) )
            b = self._generate_one_bead(self.start_position,
                                        self.num_nts,
                                        self.start_orientation)
            self.children.append(b)
            self.beads.append(b) # don't add orientation bead
            self._assign_particles_to_locations()
            return

        assert( np.linalg.norm(direction) > 0 )

172
173
174
175
176
177
178
        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
179
180
181
182
183
184
185
186
            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
187

cmaffeo2's avatar
cmaffeo2 committed
188
            b = self._generate_one_bead(pos,nts,orientation)
189
            self.children.append(b)
cmaffeo2's avatar
cmaffeo2 committed
190
191
192
193
194
            self.beads.append(b) # don't add orientation bead

            if "orientation_bead" in b.__dict__:
                o = b.orientation_bead
                self.children.append(o)
195
                self.add_bond(b,o, Segment.orientation_bond, exclude=True)
cmaffeo2's avatar
cmaffeo2 committed
196
                
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
            # 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
217
218
                 local_twist = False,
                 num_turns = None,
219
                 start_orientation = None):
cmaffeo2's avatar
cmaffeo2 committed
220
221
222
        
        self.helical_rise = 10.44
        self.distance_per_nt = 3.4
223
224
225
226
227
        Segment.__init__(self, name, num_nts, 
                         start_position,
                         end_position, 
                         segment_model)

cmaffeo2's avatar
cmaffeo2 committed
228
229
230
231
232
233
234
235
236
        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

237
238
239
240
241
242
243
244
245
        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
246
    def connect_start5(self, end3, type_="intrahelical", force_connection=False):
247
248
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
249
250
        self._connect_ends( self.start5, end3, type_, force_connection = force_connection )
    def connect_start3(self, end5, type_="intrahelical", force_connection=False):
251
252
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
253
254
        self._connect_ends( self.start3, end5, type_, force_connection = force_connection )
    def connect_end3(self, end5, type_="intrahelical", force_connection=False):
255
256
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
257
258
        self._connect_ends( self.end3, end5, type_, force_connection = force_connection )
    def connect_end5(self, end3, type_="intrahelical", force_connection=False):
259
260
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
261
        self._connect_ends( self.end5, end3, type_, force_connection = force_connection )
262
263
264

        
    ## Real work
265
    def _connect_ends(self, end1, end2, type_, force_connection):
266
267
268
269
270
271
272

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

273
        end1.container._connect( end2.container, Connection( end1, end2, type_=type_ ) )
274

cmaffeo2's avatar
cmaffeo2 committed
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
    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

293
294

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

    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
357
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
358
359
        # if (self.num_nts // max_nucleotides_per_bead) + 1 <= 1:
        #     pdb.set_trace()
cmaffeo2's avatar
cmaffeo2 committed
360
361
362
        return (self.num_nts // max_nucleotides_per_bead) + 1

    def _generate_one_bead(self, pos, nts, orientation=None):
363
364
365
366
367
368
369
370
371
372
373
374
        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
375
376
377
    def __init__(self, segments=[], local_twist=True,
                 max_basepairs_per_bead=7,
                 max_nucleotides_per_bead=4,
378
379
380
                 dimensions=(1000,1000,1000), temperature=291,
                 timestep=50e-6, cutoff=50, 
                 decompPeriod=10000, pairlistDistance=None, 
381
382
383
                 nonbondedResolution=0,DEBUG=0):
        self.DEBUG = DEBUG
        if DEBUG > 0: print("Building ARBD Model")
384
385
386
387
388
        ArbdModel.__init__(self,segments,
                           dimensions, temperature, timestep, cutoff, 
                           decompPeriod, pairlistDistance=None,
                           nonbondedResolution=0)

cmaffeo2's avatar
cmaffeo2 committed
389
390
        # self.max_basepairs_per_bead = max_basepairs_per_bead     # dsDNA
        # self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
391
392
393
394
        self.children = self.segments = segments

        self._bonded_potential = dict() # cache bonded potentials

cmaffeo2's avatar
cmaffeo2 committed
395
        self._generate_bead_model(segments, max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
396
397
398

        self.useNonbondedScheme( nbDnaScheme )

399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
    
    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 )
426
427
                    ## TODO: make this code more robust
                    ## TODO: find bug that makes things be out of order
428
                    try:
429
430
431
                        b0 = s1.get_beads_before_bead(b1,1)
                        assert(len(b0) == 1)
                        b0 = b0[0]
432
433
434
435
436
                        assert( b0 is not None )
                        ret.append( [b0,b1,b2] )
                    except:
                        ...
                    try:
437
438
439
                        b0 = s1.get_beads_after_bead(b1,1)
                        assert(len(b0) == 1)
                        b0 = b0[0]
440
441
442
443
444
                        assert( b0 is not None )
                        ret.append( [b2,b1,b0] )
                    except:
                        ...
                    try:
445
446
447
                        b3 = s2.get_beads_before_bead(b2,1)
                        assert(len(b3) == 1)
                        b3 = b3[0]
448
449
450
451
452
                        assert( b3 is not None )
                        ret.append( [b3,b2,b1] )
                    except:
                        ...
                    try:
453
454
455
                        b3 = s2.get_beads_after_bead(b2,1)
                        assert(len(b3) == 1)
                        b3 = b3[0]
456
457
458
459
460
461
                        assert( b3 is not None )
                        ret.append( [b1,b2,b3] )
                    except:
                        ...
        return ret

462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
    # 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
481
    def _get_potential(self, type_, kSpring, d, max_potential = None):
482
483
484
        key = (type_,kSpring,d)
        if key not in self._bonded_potential:
            if type_ == "bond":
485
                self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,500), max_potential=max_potential)
486
            elif type_ == "angle":
cmaffeo2's avatar
cmaffeo2 committed
487
488
                self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
                # , resolution = 1, maxForce=0.1)
489
            elif type_ == "dihedral":
cmaffeo2's avatar
cmaffeo2 committed
490
                self._bonded_potential[key] = HarmonicDihedral(kSpring,d, max_potential=max_potential)
491
492
493
494
495
496
497
            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
498
499
500
501
    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)
502
503


cmaffeo2's avatar
cmaffeo2 committed
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
    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
        

531
532
    def _generate_bead_model(self, segments,
                             max_basepairs_per_bead = 7,
cmaffeo2's avatar
cmaffeo2 committed
533
534
535
                             max_nucleotides_per_bead = 4,
                             local_twist=False
    ):
536
537
538


        """ Generate beads """
539
        if self.DEBUG: print("Generating beads")
540
        for s in segments:
cmaffeo2's avatar
cmaffeo2 committed
541
542
543
            if local_twist:
                s.local_twist = True
            s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
544
545
546
547
548
549
550
551

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

        """ Reassign bead types """
552
        if self.DEBUG: print("Assigning bead types")
553
554
555
        beadtype_s = dict()
        for segment in segments:
            for b in segment:
556
                b.num_nts = np.around( b.num_nts, decimals=1 )
557
558
559
560
561
562
563
564
565
                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
566
567
                    elif key[0] == "O":
                        t.__dict__["nts"] = b.num_nts
568
569
                    else:
                        raise Exception("TODO")
cmaffeo2's avatar
cmaffeo2 committed
570
                    # print(t.nts)
571
                    t.name = t.name + "%03d" % (10*t.nts)
572
573
574
                    beadtype_s[key] = b.type_ = t


575
        """ Add intrahelical bond potentials """
576
        if self.DEBUG: print("Adding intrahelical bond potentials")
577
        dists = dict()          # built for later use
578
579
580
        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
581
            parent = self._getParent(b1,b2)
582
583
584
585
            if b1.parent == b2.parent:
                sep = 0.5*(b1.num_nts+b2.num_nts)
            else:
                sep = 1
cmaffeo2's avatar
cmaffeo2 committed
586

587
            conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
588
            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D":
589
                elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
590
                d = 3.4*sep
591
                k = conversion*elastic_modulus/d
592
            else:
593
594
                ## TODO: get better numbers our ssDNA model
                elastic_modulus = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
595
                d = 5*sep
596
                k = conversion*elastic_modulus/d
597
                # print(sep,d,k)
598
599
600
601
602
603
604
605
              
            if b1 not in dists:
                dists[b1] = []
            if b2 not in dists:
                dists[b2] = []
            dists[b1].append([b2,sep])
            dists[b2].append([b1,sep])

606
607
608
609

            # if not (b1.type_.name[0] == "D" and b2.type_.name[0] == "D"):
            #     continue 

610
            # dists[[b1,b2]] = dists[[b2,b1]] = sep
611
612
            bond = self.get_bond_potential(k,d)
            parent.add_bond( b1, b2, bond, exclude=True )
613
614

        """ Add intrahelical angle potentials """
615
        if self.DEBUG: print("Adding intrahelical angle potentials")
616
        for b1,b2,b3 in self._get_intrahelical_angle_beads():
617
            
618
619
620
621
622
623
624
625
626
627
            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
628
            parent = self._getParent(b1,b2,b3)
629

630
631
632
633
634
635
            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
636
637
                if local_twist:
                    k *= 0.5    # halve because orientation beads have similar springs
638
639
640
641
            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
642
643
644
            # if (self.num_nts // max_nucleotides_per_bead) + 1 <= 1:
            #     pdb.set_trace()

645
646
647
648
            angle = self.get_angle_potential(k,180)
            parent.add_angle( b1, b2, b3, angle )

        """ Add intrahelical exclusions """
649
        if self.DEBUG: print("Adding intrahelical exclusions")
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
        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])] )
665
666

        if self.DEBUG: print("Adding %d exclusions" % len(exclusions))
667
        for b1,b2 in exclusions:
cmaffeo2's avatar
cmaffeo2 committed
668
            parent = self._getParent(b1,b2)
669
            parent.add_exclusion( b1, b2 )
cmaffeo2's avatar
cmaffeo2 committed
670
671
672

        """ Twist potentials """
        if local_twist:
673
            if self.DEBUG: print("Adding twist potentials")
cmaffeo2's avatar
cmaffeo2 committed
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708

            ## 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)

709
        """ Add connection potentials """
710
711
712
713
714
715
716
        for s in segments:
            for c in s.connections:
                if c.type_ == "terminal_crossover":
                    if c.A.container == s:
                        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)
717
718
719
720
721
722
723
724
725

    # 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
    #         ...