segmentmodel.py 32.2 KB
Newer Older
1
2
import numpy as np
from arbdmodel import PointParticle, ParticleType, Group, ArbdModel
3
from coords import rotationAboutAxis, quaternion_from_matrix, quaternion_to_matrix
4
5
6
7
from nonbonded import *
from copy import copy, deepcopy
from nbPot import nbDnaScheme

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

12
import types
13
14

"""
cmaffeo2's avatar
cmaffeo2 committed
15
16
17
TODO:
 - document
 - 
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
61
62


"""

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
63
64
65
66
67
68
    orientation_particle = ParticleType("O",
                                        diffusivity = 100,
                                        mass = 300,
                                        radius = 1,
                                    )
    # orientation_bond = HarmonicBond(10,2)
69
    orientation_bond = HarmonicBond(30,1.5, rRange = (0,500) )
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

    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
85
86
87
88
89
        self.start_orientation = None
        self.twist_per_nt = 0

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

90
91
92
        self._bead_model_generation = 0    # TODO: remove?
        self.segment_model = segment_model # TODO: remove?

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

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
        ## Set up interpolation for positions
        a = np.array([self.start_position,self.end_position]).T
        tck, u = interpolate.splprep( a, u=[0,1], s=0, k=1)
        self.position_spline_params = tck


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

    def contour_to_tangent(self,s):
        t = interpolate.splev( s, self.position_spline_params, der=1 )
        if len(t) > 1: t = np.array(t).T
        return t
        

    def contour_to_orientation(self,s):
        if self.start_orientation is not None:
            # axis = self.start_orientation.dot( np.array((0,0,1)) )
            if self.quaternion_spline_params is None:
                axis = self.contour_to_tangent(s)
                # orientation = rotationAboutAxis( axis, s*self.twist_per_nt*self.num_nts, normalizeAxis=False )
                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:
                q = interpolate.splev(s, self.quaternion_spline_params)
                orientation = quaternion_to_matrix(q)
        else:
            orientation = None
        return orientation


cmaffeo2's avatar
cmaffeo2 committed
134
135
136
137
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead ):
        raise NotImplementedError

    def _generate_one_bead(self, pos, nts, orientation=None):
138
139
140
141
142
143
144
        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
145
        ## Assume that consecutive beads in self.beads are bonded
146
        ret = []
cmaffeo2's avatar
cmaffeo2 committed
147
148
        for i in range(len(self.beads)-number+1):
            tmp = [self.beads[i+j] for j in range(0,number)]
149
150
151
            ret.append( tmp )
        return ret

152
    def get_beads_before_bead(self, bead, number, inclusive=False):
cmaffeo2's avatar
cmaffeo2 committed
153
154
155
        ## Assume that consecutive beads in self.beads are bonded
        i = self.beads.index(bead)
        l = len(self.beads)
156
157
158
159
160
        if i-number < 0:
            raise Exception("Not enough beads after bead")
        
        start = 1
        if inclusive: start = 0
cmaffeo2's avatar
cmaffeo2 committed
161
        return [self.beads[i-j] for j in range(start,number)]
162

163
    def get_beads_after_bead(self, bead, number, inclusive=False):
cmaffeo2's avatar
cmaffeo2 committed
164
165
166
        ## Assume that consecutive beads in self.beads are bonded
        i = self.beads.index(bead)
        l = len(self.beads)
167
168
169
170
171
        if i+number >= l:
            raise Exception("Not enough beads after bead")
        
        start = 1
        if inclusive: start = 0
cmaffeo2's avatar
cmaffeo2 committed
172
        return [self.beads[i+i] for j in range(start,number)]
173

174
175
176
177
    # def get_bead_pairs_within(self, cutoff):
    #     for b1,b2 in self.get_all_consecutive_beads(self, number)


cmaffeo2's avatar
cmaffeo2 committed
178
    def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
179
180
181
        
        """ Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """

cmaffeo2's avatar
cmaffeo2 committed
182
        ## TODO: decide whether to remove bead_model argument
183
        ##       (currently unused)
cmaffeo2's avatar
cmaffeo2 committed
184

185
186
187
        # self._bead_model_generation += 1
        # self._bead_model_max_nts_per_bead = max_nts_per_bead

cmaffeo2's avatar
cmaffeo2 committed
188
        num_beads = self._get_num_beads( max_basepairs_per_bead, max_nucleotides_per_bead )
189
        nts_per_bead = float(self.num_nts)/num_beads
cmaffeo2's avatar
cmaffeo2 committed
190
        twist_per_bead = nts_per_bead * self.twist_per_nt
191
192
193

        last = None

194
195
196
        if num_beads <= 2:
            ## not yet implemented for dsDNA
            assert( isinstance(self, SingleStrandedSegment) )
197
198
            b = self._generate_one_bead(0.5, 
                                        self.start_position,
199
200
201
202
203
204
205
206
                                        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


207
208
209
210
211
212
        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
213
            pos = self.contour_to_position(s)
cmaffeo2's avatar
cmaffeo2 committed
214
215
216
217
218
219
220
221
            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
222

223
            b = self._generate_one_bead(s,pos,nts,orientation)
224
            self.children.append(b)
cmaffeo2's avatar
cmaffeo2 committed
225
226
227
228
229
            self.beads.append(b) # don't add orientation bead

            if "orientation_bead" in b.__dict__:
                o = b.orientation_bead
                self.children.append(o)
230
                self.add_bond(b,o, Segment.orientation_bond, exclude=True)
cmaffeo2's avatar
cmaffeo2 committed
231
                
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
            # 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
252
253
                 local_twist = False,
                 num_turns = None,
254
                 start_orientation = None):
cmaffeo2's avatar
cmaffeo2 committed
255
256
257
        
        self.helical_rise = 10.44
        self.distance_per_nt = 3.4
258
259
260
261
262
        Segment.__init__(self, name, num_nts, 
                         start_position,
                         end_position, 
                         segment_model)

cmaffeo2's avatar
cmaffeo2 committed
263
264
265
266
267
268
269
270
271
        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

272
273
274
275
276
277
278
279
        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" )

280
281
282
283
284
285
286
287
288
        ## Set up interpolation for azimuthal angles 
        a = np.array([self.start_position,self.end_position]).T
        tck, u = interpolate.splprep( a, u=[0,1], s=0, k=1)
        self.position_spline_params = tck
        
        ## TODO: initialize sensible spline for orientation
        self.quaternion_spline_params = None


289
    ## Convenience methods
290
    def connect_start5(self, end3, type_="intrahelical", force_connection=False):
291
292
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
293
294
        self._connect_ends( self.start5, end3, type_, force_connection = force_connection )
    def connect_start3(self, end5, type_="intrahelical", force_connection=False):
295
296
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
297
298
        self._connect_ends( self.start3, end5, type_, force_connection = force_connection )
    def connect_end3(self, end5, type_="intrahelical", force_connection=False):
299
300
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
301
302
        self._connect_ends( self.end3, end5, type_, force_connection = force_connection )
    def connect_end5(self, end3, type_="intrahelical", force_connection=False):
303
304
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
305
        self._connect_ends( self.end5, end3, type_, force_connection = force_connection )
306
307
308

        
    ## Real work
309
    def _connect_ends(self, end1, end2, type_, force_connection):
310
311
312
313
314
315
316

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

317
        end1.container._connect( end2.container, Connection( end1, end2, type_=type_ ) )
318

cmaffeo2's avatar
cmaffeo2 committed
319
320
321
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
        return (self.num_nts // max_basepairs_per_bead) + 1

322
    def _generate_one_bead(self, contour_position, pos, nts, orientation=None):
cmaffeo2's avatar
cmaffeo2 committed
323
324
325
326
327
328
329
        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,
330
331
332
                                  num_nts=nts, parent=self, 
                                  orientation_bead=o,
                                  contour_position=contour_position )
cmaffeo2's avatar
cmaffeo2 committed
333
334
335

        else:
            bead = PointParticle( Segment.dsDNA_particle, pos, nts,
336
337
                                  num_nts=nts, parent=self,
                                  contour_position=contour_position )
cmaffeo2's avatar
cmaffeo2 committed
338
339
        return bead

340
341

    def _assign_particles_to_locations(self):
cmaffeo2's avatar
cmaffeo2 committed
342
343
        self.start3.particle =  self.start5.particle = self.beads[0]
        self.end3.particle   =  self.end5.particle   = self.beads[-1]
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403

    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
404
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead):
405
406
        # if (self.num_nts // max_nucleotides_per_bead) + 1 <= 1:
        #     pdb.set_trace()
cmaffeo2's avatar
cmaffeo2 committed
407
408
        return (self.num_nts // max_nucleotides_per_bead) + 1

409
    def _generate_one_bead(self, contour_position, pos, nts, orientation=None):
410
        return PointParticle( Segment.ssDNA_particle, pos, nts,
411
412
                              num_nts=nts, parent=self,
                              contour_position=contour_position )
413
414
415
416
417
418
419
420
421
422

    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
423
424
425
    def __init__(self, segments=[], local_twist=True,
                 max_basepairs_per_bead=7,
                 max_nucleotides_per_bead=4,
426
427
428
                 dimensions=(1000,1000,1000), temperature=291,
                 timestep=50e-6, cutoff=50, 
                 decompPeriod=10000, pairlistDistance=None, 
429
430
431
                 nonbondedResolution=0,DEBUG=0):
        self.DEBUG = DEBUG
        if DEBUG > 0: print("Building ARBD Model")
432
433
434
435
436
        ArbdModel.__init__(self,segments,
                           dimensions, temperature, timestep, cutoff, 
                           decompPeriod, pairlistDistance=None,
                           nonbondedResolution=0)

cmaffeo2's avatar
cmaffeo2 committed
437
438
        # self.max_basepairs_per_bead = max_basepairs_per_bead     # dsDNA
        # self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
439
440
441
442
        self.children = self.segments = segments

        self._bonded_potential = dict() # cache bonded potentials

443
        self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
444
445
446

        self.useNonbondedScheme( nbDnaScheme )

447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
    
    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 )
474
475
                    ## TODO: make this code more robust
                    ## TODO: find bug that makes things be out of order
476
                    try:
477
478
479
                        b0 = s1.get_beads_before_bead(b1,1)
                        assert(len(b0) == 1)
                        b0 = b0[0]
480
481
482
483
484
                        assert( b0 is not None )
                        ret.append( [b0,b1,b2] )
                    except:
                        ...
                    try:
485
486
487
                        b0 = s1.get_beads_after_bead(b1,1)
                        assert(len(b0) == 1)
                        b0 = b0[0]
488
489
490
491
492
                        assert( b0 is not None )
                        ret.append( [b2,b1,b0] )
                    except:
                        ...
                    try:
493
494
495
                        b3 = s2.get_beads_before_bead(b2,1)
                        assert(len(b3) == 1)
                        b3 = b3[0]
496
497
498
499
500
                        assert( b3 is not None )
                        ret.append( [b3,b2,b1] )
                    except:
                        ...
                    try:
501
502
503
                        b3 = s2.get_beads_after_bead(b2,1)
                        assert(len(b3) == 1)
                        b3 = b3[0]
504
505
506
507
508
509
                        assert( b3 is not None )
                        ret.append( [b1,b2,b3] )
                    except:
                        ...
        return ret

510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
    # 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
529
    def _get_potential(self, type_, kSpring, d, max_potential = None):
530
531
532
        key = (type_,kSpring,d)
        if key not in self._bonded_potential:
            if type_ == "bond":
533
                self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,500), max_potential=max_potential)
534
            elif type_ == "angle":
cmaffeo2's avatar
cmaffeo2 committed
535
536
                self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
                # , resolution = 1, maxForce=0.1)
537
            elif type_ == "dihedral":
cmaffeo2's avatar
cmaffeo2 committed
538
                self._bonded_potential[key] = HarmonicDihedral(kSpring,d, max_potential=max_potential)
539
540
541
542
543
544
545
            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
546
547
548
549
    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)
550
551


cmaffeo2's avatar
cmaffeo2 committed
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
    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
577
578
579
580
581
582
583
584
585

    # def _update_segment_positions(self, bead_coordinates):
    #     """ Set new function for each segments functions
    #     contour_to_position and contour_to_orientation """
        
    #     dsDnaHelixNeighborDist=50
    #     dsDnaAllNeighborDist=30
    #     ssDnaHelixNeighborDist=25
    #     ssDnaAllNeighborDist=25
cmaffeo2's avatar
cmaffeo2 committed
586
        
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
    #     beads = [b in s.beads for s in self.segments]
    #     positions = np.array([b.position for b in beads])
    #     neighborhood = dict()

    #     ## Assign neighborhood to each bead
    #     for b in beads:
    #         dists = b.position[np.newaxis,:] - positions
    #         dists = np.linalg.norm(dists, axis=-1)
    #         neighborhood[b] = np.where( dists < 50 )

    """ Mapping between different resolution models """
    def _clear_beads(self):
        for s in self.segments:
            s.clear_all()
            s.beads = []
        self.clear_all(keep_children=True)

    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
607

608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
        for s in self.segments:
            beads = [b for b in s.beads]
            ids = [b.idx for b in beads]
            
            """ Get positions """
            positions = bead_coordinates[ids,:].T
            contours = [b.contour_position for b in beads]
            tck, u = interpolate.splprep( positions, u=contours, s=0, )
            
            s.position_spline_params = tck

            """ Get twist """
            if 'orientation_bead' in beads[0].__dict__:
                tangents = s.contour_to_tangent(contours)
                quats = []
                for b,t in zip(beads,tangents):
                    o = b.orientation_bead
                    angleVec = o.position - b.position
                    angleVec = angleVec - angleVec.dot(t)*t
                    angleVec = angleVec/np.linalg.norm(angleVec)
                    y = np.cross(t,angleVec)
                    quats.append( quaternion_from_matrix( np.array([t,y,angleVec])) )
                quats = np.array(quats)
                tck, u = interpolate.splprep( quats.T, u=contours, s=0, )
                s.quaternion_spline_params = tck


            ## TODO: set twist

    def _generate_bead_model(self,
638
                             max_basepairs_per_bead = 7,
cmaffeo2's avatar
cmaffeo2 committed
639
                             max_nucleotides_per_bead = 4,
640
                             local_twist=False):
641

642
        segments = self.segments
643
644

        """ Generate beads """
645
        if self.DEBUG: print("Generating beads")
646
        for s in segments:
cmaffeo2's avatar
cmaffeo2 committed
647
648
649
            if local_twist:
                s.local_twist = True
            s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
650
651
652
653
654
655
656
657

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

        """ Reassign bead types """
658
        if self.DEBUG: print("Assigning bead types")
659
660
661
        beadtype_s = dict()
        for segment in segments:
            for b in segment:
662
                b.num_nts = np.around( b.num_nts, decimals=1 )
663
664
665
666
667
668
669
670
671
                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
672
673
                    elif key[0] == "O":
                        t.__dict__["nts"] = b.num_nts
674
675
                    else:
                        raise Exception("TODO")
cmaffeo2's avatar
cmaffeo2 committed
676
                    # print(t.nts)
677
                    t.name = t.name + "%03d" % (10*t.nts)
678
679
680
                    beadtype_s[key] = b.type_ = t


681
        """ Add intrahelical bond potentials """
682
        if self.DEBUG: print("Adding intrahelical bond potentials")
683
        dists = dict()          # built for later use
684
685
686
        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
687
            parent = self._getParent(b1,b2)
688
689
690
691
            if b1.parent == b2.parent:
                sep = 0.5*(b1.num_nts+b2.num_nts)
            else:
                sep = 1
cmaffeo2's avatar
cmaffeo2 committed
692

693
            conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
694
            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D":
695
                elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
696
                d = 3.4*sep
697
                k = conversion*elastic_modulus/d
698
            else:
699
700
                ## TODO: get better numbers our ssDNA model
                elastic_modulus = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
701
                d = 5*sep
702
                k = conversion*elastic_modulus/d
703
                # print(sep,d,k)
704
705
706
707
708
709
710
711
              
            if b1 not in dists:
                dists[b1] = []
            if b2 not in dists:
                dists[b2] = []
            dists[b1].append([b2,sep])
            dists[b2].append([b1,sep])

712
713
714
715

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

716
            # dists[[b1,b2]] = dists[[b2,b1]] = sep
717
718
            bond = self.get_bond_potential(k,d)
            parent.add_bond( b1, b2, bond, exclude=True )
719
720

        """ Add intrahelical angle potentials """
721
        if self.DEBUG: print("Adding intrahelical angle potentials")
722
        for b1,b2,b3 in self._get_intrahelical_angle_beads():
723
            
724
725
726
727
728
729
730
731
732
733
            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
734
            parent = self._getParent(b1,b2,b3)
735

736
737
738
739
740
741
            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
742
743
                if local_twist:
                    k *= 0.5    # halve because orientation beads have similar springs
744
745
746
747
            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
748
749
750
            # if (self.num_nts // max_nucleotides_per_bead) + 1 <= 1:
            #     pdb.set_trace()

751
752
753
754
            angle = self.get_angle_potential(k,180)
            parent.add_angle( b1, b2, b3, angle )

        """ Add intrahelical exclusions """
755
        if self.DEBUG: print("Adding intrahelical exclusions")
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
        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])] )
771
772

        if self.DEBUG: print("Adding %d exclusions" % len(exclusions))
773
        for b1,b2 in exclusions:
cmaffeo2's avatar
cmaffeo2 committed
774
            parent = self._getParent(b1,b2)
775
            parent.add_exclusion( b1, b2 )
cmaffeo2's avatar
cmaffeo2 committed
776
777
778

        """ Twist potentials """
        if local_twist:
779
            if self.DEBUG: print("Adding twist potentials")
cmaffeo2's avatar
cmaffeo2 committed
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814

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

815
        """ Add connection potentials """
816
817
818
819
820
821
822
        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)
823

824
825
826
        self._updateParticleOrder()


827
828
829
830
831
832
833
834
    # 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
    #         ...