segmentmodel.py 41.7 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
import pdb
14
"""
cmaffeo2's avatar
cmaffeo2 committed
15
TODO:
16
 - fix handling of crossovers for atomic representation
17
18
19
20
 - map to atomic representation
 - remove performance bottlenecks
 - test for large systems
 - assign sequence
21
 - document
22
23
24
25
26
"""

class Location():
    """ Site for connection within an object """
    def __init__(self, container, address, type_):
27
        ## TODO: remove cyclic references(?)
28
        self.container = container
29
        self.address = address  # represents position along contour length in segments
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
        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

cmaffeo2's avatar
cmaffeo2 committed
48
    def get_connections_and_locations(self, type_=None, exclude=[]):
49
50
51
52
        """ Returns a list with each entry of the form:
            connection, location_in_self, location_in_other """
        ret = []
        for c in self.connections:
cmaffeo2's avatar
cmaffeo2 committed
53
            if type_ is None or c.type_ == type_ and type_ not in exclude:
54
                if   c.A.container is self:
55
                    ret.append( [c, c.A, c.B] )
56
                elif c.B.container is self:
57
58
59
60
61
                    ret.append( [c, c.B, c.A] )
                else:
                    raise Exception("Object contains connection that fails to refer to object")
        return ret

62
63
64
    def _connect(self, other, connection):
        self.connections.append(connection)
        other.connections.append(connection)
65
66
    # def _find_connections(self, loc):
    #     return [c for c in self.connections if c.A == loc or c.B == loc]
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

class SegmentParticle(PointParticle):
    def __init__(self, type_, position, name="A", segname="A", **kwargs):
        self.contour_position = None
        PointParticle.__init__(self, type_, position, name=name, segname=segname, **kwargs)
        self.intrahelical_neighbors = []
        self.other_neighbors = []

    def get_contour_position(self,seg):
        assert( isinstance(seg,Segment) )
        if seg == self.parent:
            return self.contour_position
        else:
            for c,A,B in self.parent.get_connections_and_locations():
                if A.particle is self and B.container is seg:
                    nt = np.abs( (self.contour_position - A.address)*A.container.num_nts )
                    if B.address < 0.5:
                        return B.address-nt/seg.num_nts
                    else:
                        return B.address+nt/seg.num_nts
            print("")
            for c,A,B in self.parent.get_connections_and_locations():
                print("  ",c.type_)
                print(A,B)
                print(A.particle,self)
                print(B.container,seg)
            print("")
            import pdb
            pdb.set_trace()
            raise Exception("Did not find location for particle {} in Segment {}".format(self,seg))
97
98
99
100
101
102
103
104
105
106
107
108
       
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
109
110
111
112
113
    orientation_particle = ParticleType("O",
                                        diffusivity = 100,
                                        mass = 300,
                                        radius = 1,
                                    )
114

cmaffeo2's avatar
cmaffeo2 committed
115
    # orientation_bond = HarmonicBond(10,2)
116
    orientation_bond = HarmonicBond(30,1.5, rRange = (0,500) )
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

    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
132
133
134
135
136
        self.start_orientation = None
        self.twist_per_nt = 0

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

137
138
139
        self._bead_model_generation = 0    # TODO: remove?
        self.segment_model = segment_model # TODO: remove?

cmaffeo2's avatar
cmaffeo2 committed
140
        self.num_nts = int(num_nts)
141
142
143
144
145
        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

146
147
148
149
150
        ## 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

151
152
153
154
155
    def clear_all(self):
        Group.clear_all(self)  # TODO: use super?
        self.beads = []
        for c,loc,other in self.get_connections_and_locations():
            loc.particle = None
156
157
158
159
160
161
162
163

    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 )
164
165
        t = (t / np.linalg.norm(t,axis=0))
        return t.T
166
167
168
169
170
171
172
        

    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)
173
                orientation = rotationAboutAxis( axis, s*self.twist_per_nt*self.num_nts, normalizeAxis=True )
174
175
176
177
178
179
180
181
182
183
184
                ## 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
185
186
187
    def _get_num_beads(self, max_basepairs_per_bead, max_nucleotides_per_bead ):
        raise NotImplementedError

188
    def _generate_one_bead(self, contour_position, nts):
189
190
        raise NotImplementedError

191
192
193
194
195
196
197
    def get_nearest_bead(self, contour_position):
        if len(self.beads) < 1: return None
        cs = np.array([b.contour_position for b in self.beads]) # TODO: cache
        # TODO: include connected beads?
        i = np.argmin((cs - contour_position)**2)

        return self.beads[i]
198
199
200

    def get_all_consecutive_beads(self, number):
        assert(number >= 1)
cmaffeo2's avatar
cmaffeo2 committed
201
        ## Assume that consecutive beads in self.beads are bonded
202
        ret = []
cmaffeo2's avatar
cmaffeo2 committed
203
204
        for i in range(len(self.beads)-number+1):
            tmp = [self.beads[i+j] for j in range(0,number)]
205
            ret.append( tmp )
206
        return ret   
207

208
209
210
    def _add_bead(self,b,set_contour=False):
        if set_contour:
            b.contour_position = b.get_contour_position(self)
211
        
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
        # assert(b.parent is None)
        if b.parent is not None:
            b.parent.children.remove(b)
        b.parent = self
        self.children.append(b)
        self.beads.append(b) # don't add orientation bead
        if "orientation_bead" in b.__dict__: # TODO: think of a cleaner approach
            o = b.orientation_bead
            o.contour_position = b.contour_position
            if o.parent is not None:
                o.parent.children.remove(o)
            o.parent = self
            self.children.append(o)
            self.add_bond(b,o, Segment.orientation_bond, exclude=True)

    def _rebuild_children(self, new_children):
        # print("_rebuild_children on %s" % self.name)
        old_children = self.children
        old_beads = self.beads
        self.children = []
        self.beads = []
        for b in new_children:
            self.beads.append(b)
            self.children.append(b)
            if "orientation_bead" in b.__dict__: # TODO: think of a cleaner approach
                self.children.append(b.orientation_bead)
        assert(len(old_children) == len(self.children))
        assert(len(old_beads) == len(self.beads))
240

cmaffeo2's avatar
cmaffeo2 committed
241
    def _generate_beads(self, bead_model, max_basepairs_per_bead, max_nucleotides_per_bead):
242
243
        
        """ Generate beads (positions, types, etcl) and bonds, angles, dihedrals, exclusions """
cmaffeo2's avatar
cmaffeo2 committed
244
        ## TODO: decide whether to remove bead_model argument
245
        ##       (currently unused)
cmaffeo2's avatar
cmaffeo2 committed
246

247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
        ## First find points between-which beads must be generated
        locs = sorted([loc for c,loc,other_loc in self.get_connections_and_locations()], key=lambda l:l.address)
        existing_beads = [l.particle for l in locs if l.particle is not None]
        for b in existing_beads:
            assert(b.parent is not None)

        ## Add ends if they don't exist yet
        ## TODO: what about 1 nt segments?
        if len(existing_beads) == 0 or existing_beads[0].get_contour_position(self) > 0:
            b = self._generate_one_bead(0, 0)
            existing_beads = [b] + existing_beads
        if existing_beads[-1].get_contour_position(self) < 1:
            b = self._generate_one_bead(1, 0)
            existing_beads.append(b)
        assert(len(existing_beads) > 1)

        ## Walk through existing_beads, add beads between
        tmp_children = []       # build list of children in nice order
265
        last = None
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
        for I in range(len(existing_beads)-1):
            eb1,eb2 = [existing_beads[i] for i in (I,I+1)]

            # print(" %s working on %d to %d" % (self.name, eb1.position[2], eb2.position[2]))
            e_ds = eb2.get_contour_position(self) - eb1.get_contour_position(self)
            num_beads = self._get_num_beads( e_ds, max_basepairs_per_bead, max_nucleotides_per_bead )
            ds = e_ds / (num_beads+1)
            nts = ds*self.num_nts
            eb1.num_nts += 0.5*nts
            eb2.num_nts += 0.5*nts

            ## Add beads
            if eb1.parent == self:
                tmp_children.append(eb1)

            s0 = eb1.get_contour_position(self)
            if last is not None:
                last.intrahelical_neighbors.append(eb1)
                eb1.intrahelical_neighbors.append(last)
            last = eb1
            for j in range(num_beads):
                s = ds*(j+1) + s0
                b = self._generate_one_bead(s,nts)

                last.intrahelical_neighbors.append(b)
                b.intrahelical_neighbors.append(last)
                last = b
                tmp_children.append(b)

        last.intrahelical_neighbors.append(eb2)
        eb2.intrahelical_neighbors.append(last)

        if eb2.parent == self:
            tmp_children.append(eb2)
        self._rebuild_children(tmp_children)
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316

    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
317
318
                 local_twist = False,
                 num_turns = None,
319
                 start_orientation = None):
cmaffeo2's avatar
cmaffeo2 committed
320
321
322
        
        self.helical_rise = 10.44
        self.distance_per_nt = 3.4
323
324
325
326
327
        Segment.__init__(self, name, num_nts, 
                         start_position,
                         end_position, 
                         segment_model)

cmaffeo2's avatar
cmaffeo2 committed
328
329
330
331
332
333
334
335
336
        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

337
338
        self.nicks = []

339
        self.start = self.start5 = Location( self, address=0, type_= "end5" )
340
341
        self.start3 = Location( self, address=0, type_ = "end3" )

342
343
        self.end = self.end5 = Location( self, address=1, type_= "end5" )
        self.end3 = Location( self, address=1, type_ = "end3" )
344

345
346
347
348
349
350
351
352
353
        ## 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


354
    ## Convenience methods
355
    ## TODO: add errors if unrealistic connections are made
356
    def connect_start5(self, end3, type_="intrahelical", force_connection=False):
357
358
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
359
360
        self._connect_ends( self.start5, end3, type_, force_connection = force_connection )
    def connect_start3(self, end5, type_="intrahelical", force_connection=False):
361
362
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
363
364
        self._connect_ends( self.start3, end5, type_, force_connection = force_connection )
    def connect_end3(self, end5, type_="intrahelical", force_connection=False):
365
366
        if isinstance(end5, SingleStrandedSegment):
            end5 = end5.end5
367
368
        self._connect_ends( self.end3, end5, type_, force_connection = force_connection )
    def connect_end5(self, end3, type_="intrahelical", force_connection=False):
369
370
        if isinstance(end3, SingleStrandedSegment):
            end3 = end3.end3
371
        self._connect_ends( self.end5, end3, type_, force_connection = force_connection )
372

cmaffeo2's avatar
cmaffeo2 committed
373
    def crossover(self, nt, other, other_nt, strands_fwd=[True,False]):
cmaffeo2's avatar
cmaffeo2 committed
374
375
376
377
378
379
380
381
        """ Add a crossover between two helices """
        ## Validate other, nt, other_nt
        ##   TODO

        ## Create locations, connections and add to segments
        c = nt/self.num_nts
        print(self.name,nt,c)
        assert(c >= 0 and c <= 1)
cmaffeo2's avatar
cmaffeo2 committed
382
        loc = Location( self, address=c, type_ = strands_fwd[0] )
cmaffeo2's avatar
cmaffeo2 committed
383
384
385
386

        c = other_nt/other.num_nts
        print(other.name,other_nt,c)
        assert(c >= 0 and c <= 1)
cmaffeo2's avatar
cmaffeo2 committed
387
388
        other_loc = Location( other, address=c, type_ = strands_fwd[1] )
        self._connect(other, Connection( loc, other_loc, type_="crossover" ))
cmaffeo2's avatar
cmaffeo2 committed
389

390
    ## Real work
391
    def _connect_ends(self, end1, end2, type_, force_connection):
392
        ## TODO remove self?
393
394
395
396
397
        ## validate the input
        for end in (end1, end2):
            assert( isinstance(end, Location) )
            assert( end.type_ in ("end3","end5") )
        assert( end1.type_ != end2.type_ )
398
        ## Create and add connection
399
        end1.container._connect( end2.container, Connection( end1, end2, type_=type_ ) )
400

401
402
    def _get_num_beads(self, contour, max_basepairs_per_bead, max_nucleotides_per_bead):
        return int(contour*self.num_nts // max_basepairs_per_bead)
cmaffeo2's avatar
cmaffeo2 committed
403

404
405
    def _generate_one_bead(self, contour_position, nts):
        pos = self.contour_to_position(contour_position)
cmaffeo2's avatar
cmaffeo2 committed
406
        if self.local_twist:
407
            orientation = self.contour_to_orientation(contour_position)
cmaffeo2's avatar
cmaffeo2 committed
408
            opos = pos + orientation.dot( np.array((Segment.orientation_bond.r0,0,0)) )
409
410
411
412
413
414
            o = SegmentParticle( Segment.orientation_particle, opos, nts,
                                 num_nts=nts, parent=self )
            bead = SegmentParticle( Segment.dsDNA_particle, pos, nts,
                                    num_nts=nts, parent=self, 
                                    orientation_bead=o,
                                    contour_position=contour_position )
cmaffeo2's avatar
cmaffeo2 committed
415
416

        else:
417
418
419
420
            bead = SegmentParticle( Segment.dsDNA_particle, pos, nts,
                                    num_nts=nts, parent=self,
                                    contour_position=contour_position )
        self._add_bead(bead)
cmaffeo2's avatar
cmaffeo2 committed
421
422
        return bead

423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463

    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" )
464
        self.end = self.end3 = Location( self, address=1, type_ = "end3" )
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482

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

483
484
    def _get_num_beads(self, contour, max_basepairs_per_bead, max_nucleotides_per_bead):
        return int(contour*self.num_nts // max_nucleotides_per_bead)
cmaffeo2's avatar
cmaffeo2 committed
485

486
487
    def _generate_one_bead(self, contour_position, nts):
        pos = self.contour_to_position(contour_position)
488
489
490
491
492
        b = SegmentParticle( Segment.ssDNA_particle, pos, nts,
                             num_nts=nts, parent=self,
                             contour_position=contour_position )
        self._add_bead(b)
        return b
493
494
495
496
497
498

    def _generate_atomic(self, atomic_model):
        ...
    

class SegmentModel(ArbdModel):
cmaffeo2's avatar
cmaffeo2 committed
499
500
501
    def __init__(self, segments=[], local_twist=True,
                 max_basepairs_per_bead=7,
                 max_nucleotides_per_bead=4,
502
503
504
                 dimensions=(1000,1000,1000), temperature=291,
                 timestep=50e-6, cutoff=50, 
                 decompPeriod=10000, pairlistDistance=None, 
505
506
507
                 nonbondedResolution=0,DEBUG=0):
        self.DEBUG = DEBUG
        if DEBUG > 0: print("Building ARBD Model")
508
509
510
511
512
        ArbdModel.__init__(self,segments,
                           dimensions, temperature, timestep, cutoff, 
                           decompPeriod, pairlistDistance=None,
                           nonbondedResolution=0)

cmaffeo2's avatar
cmaffeo2 committed
513
514
        # self.max_basepairs_per_bead = max_basepairs_per_bead     # dsDNA
        # self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
515
516
517
518
        self.children = self.segments = segments

        self._bonded_potential = dict() # cache bonded potentials

519
        self._generate_bead_model( max_basepairs_per_bead, max_nucleotides_per_bead, local_twist)
520
521
522

        self.useNonbondedScheme( nbDnaScheme )

523

cmaffeo2's avatar
cmaffeo2 committed
524
    def get_connections(self,type_=None,exclude=[]):
525
526
527
528
        """ Find all connections in model, without double-counting """
        added=set()
        ret=[]
        for s in self.segments:
cmaffeo2's avatar
cmaffeo2 committed
529
            items = [e for e in s.get_connections_and_locations(type_,exclude=exclude) if e[0] not in added]
530
531
532
            added.update([e[0] for e in items])
            ret.extend( items )
        return ret
533
    
534
    def _recursively_get_beads_within_bonds(self,b1,bonds,done=[]):
535
        ret = []
536
537
538
539
540
541
542
543
544
        done = list(done)
        done.append(b1)
        if bonds == 0:
            return [[]]

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

547
548
549
550
551
    def _get_intrahelical_beads(self,num=2):
        ## TODO: add check that this is not called before adding intrahelical_neighbors in _generate_bead_model

        assert(num >= 2)

552
553
        ret = []
        for s in self.segments:
554
555
556
557
558
            for b1 in s.beads:
                for bead_list in self._recursively_get_beads_within_bonds(b1, num-1):
                    assert(len(bead_list) == num-1)
                    if b1.idx < bead_list[-1].idx: # avoid double-counting
                        ret.append([b1]+bead_list)
559
560
        return ret

561

562
563
    def _get_intrahelical_angle_beads(self):
        return self._get_intrahelical_beads(num=3)
564

cmaffeo2's avatar
cmaffeo2 committed
565
    def _get_potential(self, type_, kSpring, d, max_potential = None):
566
567
568
        key = (type_,kSpring,d)
        if key not in self._bonded_potential:
            if type_ == "bond":
569
                self._bonded_potential[key] = HarmonicBond(kSpring,d, rRange=(0,500), max_potential=max_potential)
570
            elif type_ == "angle":
cmaffeo2's avatar
cmaffeo2 committed
571
572
                self._bonded_potential[key] = HarmonicAngle(kSpring,d, max_potential=max_potential)
                # , resolution = 1, maxForce=0.1)
573
            elif type_ == "dihedral":
cmaffeo2's avatar
cmaffeo2 committed
574
                self._bonded_potential[key] = HarmonicDihedral(kSpring,d, max_potential=max_potential)
575
576
577
578
579
580
581
            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
582
583
584
585
    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)
586
587


cmaffeo2's avatar
cmaffeo2 committed
588
589
590
    def _getParent(self, *beads ):
        ## TODO: test
        if np.all( [b1.parent == b2.parent 
591
                    for b1,b2 in zip(beads[:-1],beads[1:])] ):
cmaffeo2's avatar
cmaffeo2 committed
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
            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
613
614
615
616
617
618
619
620
621

    # 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
622
        
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
    #     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()
        self.clear_all(keep_children=True)
638
639
640
641
642
643
644
645
646
647
        assert( len([b for b in self]) == 0 )
        locParticles = []
        # for c,A,B in self.get_connections():
        for s in self.segments:
            for c,A,B in s.get_connections_and_locations():
                for l in (A,B):
                    if l.particle is not None:
                        locParticles.append(A.particle)
        assert( len(locParticles) == 0 )
        assert( len([b for s in self.segments for b in s.beads]) == 0 )
648
649
650
651

    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
652

653
        for s in self.segments:
654
655
656
657
658
659
660
661
662
            ## TODO: may need to make spline on helices continuously differentiable
            # beads = s.beads # [b for b in s.beads]
            beads = list(set(s.beads + [A.particle for c,A,B in s.get_connections_and_locations("intrahelical")]))
            contours = [b.get_contour_position(s) for b in beads]

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

663
664
665
666
            ids = [b.idx for b in beads]
            
            """ Get positions """
            positions = bead_coordinates[ids,:].T
667
668
669
670
671
672
673

            ## TODO: determine whether having beads is an issue 
            try:
                tck, u = interpolate.splprep( positions, u=contours, s=0, k=3 )
            except:
                tck, u = interpolate.splprep( positions, u=contours, s=0, k=1 )

674
675
676
            s.position_spline_params = tck

            """ Get twist """
677
678
679
680
681
682
            cb = [e for e in cb if 'orientation_bead' in e[1].__dict__]
            beads = [b for c,b in cb] 
            contours = [c for c,b in cb] 
            ids = [b.idx for b in beads]
            # if 'orientation_bead' in beads[0].__dict__:
            if len(beads) > 3:
683
684
685
686
687
688
689
690
691
692
                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)
693
694
695
696
                try:
                    tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=3 )
                except:
                    tck, u = interpolate.splprep( quats.T, u=contours, s=0, k=1 )
697
698
699
700
701
702
                s.quaternion_spline_params = tck


            ## TODO: set twist

    def _generate_bead_model(self,
703
                             max_basepairs_per_bead = 7,
cmaffeo2's avatar
cmaffeo2 committed
704
                             max_nucleotides_per_bead = 4,
705
                             local_twist=False):
706

707
        segments = self.segments
708
709
710
        for s in segments:
            s.local_twist = local_twist
            
711
        """ Generate beads at intrahelical junctions """
712
713
714
715
716
717
718
719
720
721
722
723
724
725
        if self.DEBUG: print( "Adding intrahelical beads at junctions" )
        ## Loop through all connections, generate beads at appropriate locations
        for c,A,B in self.get_connections("intrahelical"):
            s1,s2 = [l.container for l in (A,B)]
            if isinstance(s1,DoubleStrandedSegment) and isinstance(s2,DoubleStrandedSegment):
                assert( A.particle is None )
                assert( B.particle is None )

                ## TODO: offload the work here to s1
                a1,a2 = [l.address   for l in (A,B)]
                a1,a2 = [a - (0.5/s.num_nts) if a == 0 else a + (0.5/s.num_nts) for a,s in zip((a1,a2),(s1,s2))]

                b = s1._generate_one_bead(a1,0)
                A.particle = B.particle = b
726

727
            else:
cmaffeo2's avatar
cmaffeo2 committed
728
                ## TODO fix this for ssDNA vs dsDNA (maybe a1/a2 should be calculated differently)
729
730
                a1,a2 = [l.address   for l in (A,B)]
                a1,a2 = [a - (0.5/s.num_nts) if a == 0 else a + (0.5/s.num_nts) for a,s in zip((a1,a2),(s1,s2))]
cmaffeo2's avatar
cmaffeo2 committed
731
732
733
734
                if isinstance(s2,DoubleStrandedSegment):
                    b = s2._generate_one_bead(a2,0)
                else:
                    b = s1._generate_one_bead(a1,0)
735
736
                A.particle = B.particle = b

737
        """ Generate beads at other junctions """
cmaffeo2's avatar
cmaffeo2 committed
738
739
740
741
742
743
744
        for c,A,B in self.get_connections(exclude="intrahelical"):
            s1,s2 = [l.container for l in (A,B)]
            if A.particle is not None and B.particle is not None:
                continue
            assert( A.particle is None )
            assert( B.particle is None )

745
            ## TODO: offload the work here to s1/s2 (?)
cmaffeo2's avatar
cmaffeo2 committed
746
            a1,a2 = [l.address   for l in (A,B)]
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762

            b = s1.get_nearest_bead(a1)
            if b is not None and np.abs(b.contour_position-a1)*s1.num_nts < 1:
                ## combine beads
                b.contour_position = 0.5*(b.contour_position + a1) # avg position
                A.particle = b
            else:
                A.particle = s1._generate_one_bead(a1,0)

            b = s2.get_nearest_bead(a2)
            if b is not None and np.abs(b.contour_position-a2)*s2.num_nts < 1:
                ## combine beads
                b.contour_position = 0.5*(b.contour_position + a2) # avg position
                B.particle = b
            else:
                B.particle = s2._generate_one_bead(a2,0)
cmaffeo2's avatar
cmaffeo2 committed
763

764
765
766
767
768
769
770
        """ Some tests """
        for c,A,B in self.get_connections("intrahelical"):
            for l in (A,B):
                if l.particle is None: continue
                assert( l.particle.parent is not None )

        """ Generate beads in between """
771
        if self.DEBUG: print("Generating beads")
772
        for s in segments:
cmaffeo2's avatar
cmaffeo2 committed
773
            s._generate_beads( self, max_basepairs_per_bead, max_nucleotides_per_bead )
774

775
776
777
778
779
780
781
782
783
784
785
786
787
788
        # """ Combine beads at junctions as needed """
        # for c,A,B in self.get_connections():
        #    ...

        """ Add intrahelical neighbors at connections """
        for c,A,B in self.get_connections("intrahelical"):
            b1,b2 = [l.particle for l in (A,B)]
            if b1 is b2:
                ## already handled by Segment._generate_beads
                continue
            else:
                for b in (b1,b2): assert( b is not None )
                b1.intrahelical_neighbors.append(b2)
                b2.intrahelical_neighbors.append(b1)
789
790

        """ Reassign bead types """
791
        if self.DEBUG: print("Assigning bead types")
792
793
794
        beadtype_s = dict()
        for segment in segments:
            for b in segment:
795
                b.num_nts = np.around( b.num_nts, decimals=1 )
796
797
798
799
800
801
802
803
804
                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
805
806
                    elif key[0] == "O":
                        t.__dict__["nts"] = b.num_nts
807
808
                    else:
                        raise Exception("TODO")
cmaffeo2's avatar
cmaffeo2 committed
809
                    # print(t.nts)
810
                    t.name = t.name + "%03d" % (10*t.nts)
811
812
                    beadtype_s[key] = b.type_ = t

813
814
815
        """ Update bead indices """
        self._countParticleTypes() # probably not needed here
        self._updateParticleOrder()
816

817
        """ Add intrahelical bond potentials """
818
        if self.DEBUG: print("Adding intrahelical bond potentials")
cmaffeo2's avatar
cmaffeo2 committed
819
        dists = dict()          # intrahelical distances built for later use
820
821
822
        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
823
            parent = self._getParent(b1,b2)
824
825
826

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

828
            conversion = 0.014393265 # units "pN/AA" kcal_mol/AA^2
829
            if b1.type_.name[0] == "D" and b2.type_.name[0] == "D":
830
                elastic_modulus = 1000 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
831
                d = 3.4*sep
832
                k = conversion*elastic_modulus/d
833
            else:
834
835
                ## TODO: get better numbers our ssDNA model
                elastic_modulus = 800 # pN http://markolab.bmbcb.northwestern.edu/marko/Cocco.CRP.02.pdf
836
                d = 5*sep
837
                k = conversion*elastic_modulus/d
838
                # print(sep,d,k)
839
840
              
            if b1 not in dists:
cmaffeo2's avatar
cmaffeo2 committed
841
                dists[b1] = dict()
842
            if b2 not in dists:
cmaffeo2's avatar
cmaffeo2 committed
843
844
845
846
847
                dists[b2] = dict()
            # dists[b1].append([b2,sep])
            # dists[b2].append([b1,sep])
            dists[b1][b2] = sep
            dists[b2][b1] = sep
848
849

            # dists[[b1,b2]] = dists[[b2,b1]] = sep
850
851
            bond = self.get_bond_potential(k,d)
            parent.add_bond( b1, b2, bond, exclude=True )
852
853

        """ Add intrahelical angle potentials """
854
        if self.DEBUG: print("Adding intrahelical angle potentials")
855
        for b1,b2,b3 in self._get_intrahelical_angle_beads():
856
857
            ## TODO: could be slightly smarter about sep
            sep = 0.5*b1.num_nts+b2.num_nts+0.5*b3.num_nts
cmaffeo2's avatar
cmaffeo2 committed
858
            parent = self._getParent(b1,b2,b3)
859

860
861
862
863
864
865
            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
866
867
                if local_twist:
                    k *= 0.5    # halve because orientation beads have similar springs
868
869
            else:
                ## TODO: get correct number from ssDNA model
cmaffeo2's avatar
cmaffeo2 committed
870
                k = 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/3))) * 0.00030461742; # kcal_mol/degree^2
871

872
873
874
875
            angle = self.get_angle_potential(k,180)
            parent.add_angle( b1, b2, b3, angle )

        """ Add intrahelical exclusions """
876
        if self.DEBUG: print("Adding intrahelical exclusions")
877
878
879
        beads = dists.keys()
        def _recursively_get_beads_within(b1,d,done=[]):
            ret = []
cmaffeo2's avatar
cmaffeo2 committed
880
            for b2,sep in dists[b1].items():
881
882
883
884
885
886
887
                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
888

889
890
891
        exclusions = set()
        for b1 in beads:
            exclusions.update( [(b1,b) for b in _recursively_get_beads_within(b1, 20, done=[b1])] )
892
893

        if self.DEBUG: print("Adding %d exclusions" % len(exclusions))
894
        for b1,b2 in exclusions:
cmaffeo2's avatar
cmaffeo2 committed
895
            parent = self._getParent(b1,b2)
896
            parent.add_exclusion( b1, b2 )
cmaffeo2's avatar
cmaffeo2 committed
897
898
899

        """ Twist potentials """
        if local_twist:
900
            if self.DEBUG: print("Adding twist potentials")
cmaffeo2's avatar
cmaffeo2 committed
901
902
903
904
905
906
907
908
909
910

            ## 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
cmaffeo2's avatar
cmaffeo2 committed
911
                for b2,sep in dists[b1].items():
cmaffeo2's avatar
cmaffeo2 committed
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
                    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)

936
        """ Add connection potentials """
937
        for c,A,B in self.get_connections("terminal_crossover"):
cmaffeo2's avatar
cmaffeo2 committed
938
            ## TODO: use a better description here
939
940
941
            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)
942

cmaffeo2's avatar
cmaffeo2 committed
943
        for c,A,B in self.get_connections("crossover"):
cmaffeo2's avatar
cmaffeo2 committed
944
            ## TODO: avoid double-counting for double-crossovers
cmaffeo2's avatar
cmaffeo2 committed
945
946
947
948
949
950
951
952
953
954
955
            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)

            if not local_twist:
                ## TODO
                ...
            else:
                k = (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(1)/147))) * 0.00030461742; # kcal_mol/degree^2
                if 'orientation_bead' in b1.__dict__:
                    t0 = 90 + 60
cmaffeo2's avatar
cmaffeo2 committed
956
                    if A.type_: t0 -= 120
cmaffeo2's avatar
cmaffeo2 committed
957
958
959
960
961
                    o = b1.orientation_bead
                    pot = self.get_angle_potential(k,t0)
                    self.add_angle( o,b1,b2, pot )
                else:
                    t0 = 90 + 60
cmaffeo2's avatar
cmaffeo2 committed
962
                    if B.type_: t0 -= 120
cmaffeo2's avatar
cmaffeo2 committed
963
964
965
966
967
                    o = b2.orientation_bead
                    pot = self.get_angle_potential(k,t0)
                    self.add_angle( b1,b2,o, pot )


cmaffeo2's avatar
cmaffeo2 committed
968
969
                ## TODO: fix the following to ensure that neighbors are indeed above and below
                ### preferably with a function call to b
cmaffeo2's avatar
cmaffeo2 committed
970
971
972
                u1,u2 = [b.intrahelical_neighbors[-1] for b in (b1,b2)] # TODO: fix this
                d1,d2 = [b.intrahelical_neighbors[0] for b in (b1,b2)] # TODO: fix this

cmaffeo2's avatar
cmaffeo2 committed
973
                k_fn = lambda sep: (1.0/2) * 1.5 * kT * (1.0 / (1-np.exp(-float(sep)/147))) * 0.00030461742; # kcal_mol/degree^2
cmaffeo2's avatar
cmaffeo2 committed
974
975
976
977
                t0 = 90
                if 'orientation_bead' in b1.__dict__:
                    o1 = b1.orientation_bead
                    if u2 is not None:
cmaffeo2's avatar
cmaffeo2 committed
978
979
                        k = k_fn( dists[b2][u2] )
                        pot = self.get_dihedral_potential(k,t0)
cmaffeo2's avatar
cmaffeo2 committed
980
981
                        self.add_dihedral( o1,b1,b2,u2, pot )
                    elif d2 is not None:
cmaffeo2's avatar
cmaffeo2 committed
982
983
984
                        k = k_fn( dists[b2][d2] )
                        pot = self.get_dihedral_potential(k,t0)
                        self.add_dihedral( o1,b1,b2,d2, pot )
cmaffeo2's avatar
cmaffeo2 committed
985
986
987
                if 'orientation_bead' in b2.__dict__:
                    o2 = b2.orientation_bead
                    if u1 is not None:
cmaffeo2's avatar
cmaffeo2 committed
988
989
                        k = k_fn( dists[b1][u1] )
                        pot = self.get_dihedral_potential(k,t0)
cmaffeo2's avatar
cmaffeo2 committed
990
991
                        self.add_dihedral( o2,b2,b1,u1, pot )
                    elif d1 is not None:
cmaffeo2's avatar
cmaffeo2 committed
992
993
                        k = k_fn( dists[b1][d1] )
                        pot = self.get_dihedral_potential(k,t0)
cmaffeo2's avatar
cmaffeo2 committed
994
995
996
997
998
                        self.add_dihedral( o2,b2,b1,d1, pot )

                if 'orientation_bead' in b1.__dict__ and 'orientation_bead' in b2.__dict__:
                    if u1 is not None and u2 is not None:
                        t0 = 0
cmaffeo2's avatar
cmaffeo2 committed
999
                        k = k_fn( 0.5*(dists[b1][u1]+dists[b2][u2]) )
cmaffeo2's avatar
cmaffeo2 committed
1000
                        pot = self.get_dihedral_potential(k,t0)
For faster browsing, not all history is shown. View entire blame