segmentmodel.py 165 KB
Newer Older
cmaffeo2's avatar
cmaffeo2 committed
3001
3002
            nt1,nt2 = [l.get_nt_pos() for l in (A,B)]
            is_end1, is_end2 = [nt in (0,l.container.num_nt-1) for nt,l in zip((nt1,nt2),(A,B))]
3003
3004
3005
3006
3007
3008
3009
3010
3011
            is_T_junction = (is_end1 and not is_end2) or (is_end2 and not is_end1)

            if (not is_end1) and (not is_end2):
                ## TODO?: Only apply this potential if not local_twist
                add_parallel_crossover_potential(b1,b2)

            # dotProduct = b1.parent.contour_to_tangent(b1.contour_position).dot(
            #     b2.parent.contour_to_tangent(b2.contour_position) )

3012

3013
            if local_twist:
3014
3015
                if is_T_junction:
                    """ Special case: one helix extends away from another in T-shaped junction """
3016
3017
3018
3019
3020
                    """ Add bond potential """
                    k = 1.0 * count_crossovers((b1,b2))
                    pot = self.get_bond_potential(k,12.5)
                    self.add_bond(b1,b2, pot)

3021
                    if is_end1:
cmaffeo2's avatar
cmaffeo2 committed
3022
3023
                        b1_forward = A.on_fwd_strand if nt1 == 0 else not A.on_fwd_strand
                        add_local_tee_orientation_potential(b1,b2, b1_forward, B.on_fwd_strand)
3024
                    else:
cmaffeo2's avatar
cmaffeo2 committed
3025
                        add_local_crossover_strand_orientation_potential(b1,b2, A.on_fwd_strand)
cmaffeo2's avatar
cmaffeo2 committed
3026

3027
                    if is_end2:
cmaffeo2's avatar
cmaffeo2 committed
3028
3029
                        b2_forward = B.on_fwd_strand if nt2 == 0 else not B.on_fwd_strand
                        add_local_tee_orientation_potential(b2,b1, b2_forward, A.on_fwd_strand)
3030
                    else:
cmaffeo2's avatar
cmaffeo2 committed
3031
                        add_local_crossover_strand_orientation_potential(b2,b1, B.on_fwd_strand)
3032
3033

                else:
3034
3035
3036
3037
3038
                    """ Add bond potential """
                    k = 1.0 * count_crossovers((b1,b2))
                    pot = self.get_bond_potential(k,18.5)
                    self.add_bond(b1,b2, pot)

3039
                    """ Normal case: add orientation potential """
cmaffeo2's avatar
cmaffeo2 committed
3040
3041
                    add_local_crossover_strand_orientation_potential(b1,b2, A.on_fwd_strand)
                    add_local_crossover_strand_orientation_potential(b2,b1, B.on_fwd_strand)
3042
3043
3044
3045
3046
3047
            else:
                """ Add bond potential """
                k = 1.0 * count_crossovers((b1,b2))
                pot = self.get_bond_potential(k,18.5)
                self.add_bond(b1,b2, pot)

3048
3049
3050
3051
3052

        """ Add connection potentials """
        processed_crossovers = set()
        # pdb.set_trace()
        for c,A,B in self.get_connections("sscrossover"):
cmaffeo2's avatar
cmaffeo2 committed
3053
            p1,p2 = [loc.container for loc in (A,B)]
3054
3055
3056
3057
3058

            assert(any([isinstance(p,SingleStrandedSegment) for p in (p1,p2)]))
            add_ss_crossover_potentials(c,A,B)

        for c,A,B in self.get_connections("intrahelical"):
cmaffeo2's avatar
cmaffeo2 committed
3059
            ps = [loc.container for loc in (A,B)]
3060
3061
3062
3063
3064
3065

            if any([isinstance(p,SingleStrandedSegment) for p in ps]) and \
               any([isinstance(p,DoubleStrandedSegment) for p in ps]):
                add_ss_crossover_potentials(c,A,B, add_bond=False)

        for c,A,B in sum([self.get_connections(term) for term in ("crossover","terminal_crossover")],[]):
cmaffeo2's avatar
cmaffeo2 committed
3066
            p1,p2 = [loc.container for loc in (A,B)]
3067
3068
3069
3070
3071
3072
3073

            if any([isinstance(p,SingleStrandedSegment) for p in (p1,p2)]):
                add_ss_crossover_potentials(c,A,B)
            else:
                add_crossover_potentials(c,A,B)

        ## todotodo check that this works
3074
        for crossovers in self.get_consecutive_crossovers():
cmaffeo2's avatar
cmaffeo2 committed
3075
            if local_twist: break
3076
            ## filter crossovers
cmaffeo2's avatar
cmaffeo2 committed
3077
            for i in range(len(crossovers)-2):
3078
                c1,A1,B1,dir1 = crossovers[i]
cmaffeo2's avatar
cmaffeo2 committed
3079
                c2,A2,B2,dir2 = crossovers[i+1]
3080
                s1,s2 = [l.container for l in (A1,A2)]
3081
                sep = A1.particle.get_nt_position(s1,near_address=A1.address) - A2.particle.get_nt_position(s2,near_address=A2.address)
3082
3083
                sep = np.abs(sep)

3084
3085
                assert(sep >= 0)

3086
3087
                n1,n2,n3,n4 = (B1.particle, A1.particle, A2.particle, B2.particle)

3088
3089
3090
                """
                <cos(q)> = exp(-s/Lp) = integrate( cos[x] exp(-A x^2), {x, 0, pi} ) / integrate( exp(-A x^2), {x, 0, pi} )
                """
3091
3092
3093

                ## From http://www.annualreviews.org/doi/pdf/10.1146/annurev.bb.17.060188.001405
                ##   units "3e-19 erg cm/ 295 k K" "nm" =~ 73
cmaffeo2's avatar
cmaffeo2 committed
3094
                Lp = s1.twist_persistence_length/0.34  # set semi-arbitrarily as there is a large spread in literature
3095

3096
                def get_spring(sep):
3097
                    kT = self.temperature * 0.0019872065 # kcal/mol
3098
3099
3100
3101
                    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

3102
3103
3104
                k = get_spring( max(sep,1) )
                k = (1/k + 2/k_xover_angle(sep=1))**-1

3105
                t0 = sep*s1.twist_per_nt # TODO weighted avg between s1 and s2
3106

3107
                # pdb.set_trace()
cmaffeo2's avatar
cmaffeo2 committed
3108
                if A1.on_fwd_strand: t0 -= 120
3109
                if dir1 != dir2:
cmaffeo2's avatar
cmaffeo2 committed
3110
                    A2_on_fwd = not A2.on_fwd_strand
3111
                else:
cmaffeo2's avatar
cmaffeo2 committed
3112
3113
                    A2_on_fwd = A2.on_fwd_strand
                if A2_on_fwd: t0 += 120
3114
3115
3116
3117
3118
   
                # t0 = (t0 % 360

                # if n2.idx == 0:
                #     print( n1.idx,n2.idx,n3.idx,n4.idx,k,t0,sep )
3119
                if sep == 0:
3120
3121
3122
                    # pot = self.get_angle_potential(k,t0)
                    # self.add_angle( n1,n2,n4, pot )
                    pass
3123
3124
                elif len(set((n1,n2,n3,n4))) != 4:
                    print("WARNING: skipping crossover dihedral angle potential because beads are too close")
3125
3126
3127
                else:
                    pot = self.get_dihedral_potential(k,t0)
                    self.add_dihedral( n1,n2,n3,n4, pot )
3128

3129
3130
3131
3132
        for callback in self._generate_bead_callbacks:
            callback(self)


3133
3134
        # ## remove duplicate potentials; ## TODO ensure that they aren't added twice in the first place? 
        # self.remove_duplicate_terms()
3135

3136
    def walk_through_helices(segment, direction=1, processed_segments=None):
3137
3138
3139
3140
3141
        """
    
        First and last segment should be same for circular helices
        """

3142
        assert( direction in (1,-1) )
3143
        if processed_segments == None:
3144
            processed_segments = set()
3145
3146

        def segment_is_new_helix(s):
3147
            return isinstance(s,DoubleStrandedSegment) and s not in processed_segments
3148

3149
3150
3151
        new_s = None
        s = segment
        ## iterate intrahelically connected dsDNA segments
3152
3153
        while segment_is_new_helix(s):
            conn_locs = s.get_contour_sorted_connections_and_locations("intrahelical")[::direction]
3154
            processed_segments.add(new_s)        
3155
3156
3157
3158
3159
            new_s = None
            new_dir = None
            for i in range(len(conn_locs)):
                c,A,B = conn_locs[i]
                ## TODO: handle change of direction
cmaffeo2's avatar
cmaffeo2 committed
3160
                # TODOTODO
3161
3162
3163
3164
3165
3166
3167
3168
3169
3170
                address = 1*(direction==-1)
                if A.address == address and segment_is_new_helix(B.container):
                    new_s = B.container
                    assert(B.address in (0,1))
                    new_dir = 2*(B.address == 0) - 1
                    break

            yield s,direction
            s = new_s   # will break if None
            direction = new_dir
3171
            
3172
3173
3174
3175
3176
3177
3178
3179
3180
        #     if new_s is None:
        #         break
        #     else:
        #         s = new_s
        # yield s
        ## return s


    def get_consecutive_crossovers(self):
3181
        ## TODOTODO TEST
3182
3183
3184
3185
3186
3187
3188
3189
3190
3191
3192
3193
        crossovers = []
        processed_segments = set()
        for s1 in self.segments:
            if not isinstance(s1,DoubleStrandedSegment):
                continue
            if s1 in processed_segments: continue

            s0,d0 = list(SegmentModel.walk_through_helices(s1,direction=-1))[-1]

            # s,direction = get_start_of_helices()
            tmp = []
            for s,d in SegmentModel.walk_through_helices(s0,-d0):
3194
3195
                if s == s0 and len(tmp) > 0:
                    ## end of circular helix, only add first crossover
3196
3197
3198
                    cl_list = s.get_contour_sorted_connections_and_locations("crossover")
                    if len(cl_list) > 0:
                        tmp.append( cl_list[::d][0] + [d] )
3199
3200
                else:
                    tmp.extend( [cl + [d] for cl in s.get_contour_sorted_connections_and_locations("crossover")[::d]] )
3201
3202
                processed_segments.add(s)
            crossovers.append(tmp)
3203
        return crossovers
3204

3205
3206
3207
3208
3209
3210
3211
3212
3213
3214
3215
    def set_sequence(self, sequence, force=True):
        if force:
            self.strands[0].set_sequence(sequence)
        else:
            try:
                self.strands[0].set_sequence(sequence)
            except:
                ...
        for s in self.segments:
            s.randomize_unset_sequence()

3216
3217
3218
3219
3220
3221
3222
3223
3224
3225
3226
3227
3228
3229
3230
3231
3232
3233
3234
3235
3236
3237
3238
3239
3240
3241
3242
3243
3244
    def _set_cadnano2_sequence(self, cadnano_sequence_csv_file, replace_unset=None):
        try:
            self.segments[0]._cadnano_helix
        except:
            raise Exception("Cannot set sequence of a non-cadnano model using a cadnano sequence file")

        starts = dict()
        for strand in self.strands[1:]:
            if strand.is_circular: continue
            s = strand.strand_segments[0]
            hid = s.segment._cadnano_helix
            idx = s.segment._cadnano_bp_to_zidx[int(s.start)]
            starts["{:d}[{:d}]".format(hid,idx)] = strand

        base_set = tuple('ATCG')
        with open(cadnano_sequence_csv_file) as fh:
            reader = csv.reader(fh)
            for values in reader:
                if values[0] == "Start": continue
                start,end,seq = values[:3]
                strand = starts[start]
                if replace_unset in base_set:
                    seq = [replace_unset if s == '?' else s for s in seq]
                else:
                    seq = [random.choice(base_set) if s == '?' else s for s in seq]
                if len([s for s in seq if s not in base_set]) > 0:
                    print(seq)
                strand.set_sequence(seq)

3245
    def _generate_strands(self):
3246
3247
3248
        ## clear strands
        try:
            for s in self.strands:
3249
3250
3251
3252
3253
3254
3255
3256
                try:
                    self.children.remove(s)
                except:
                    pass
        except:
            pass

        try:
3257
            for seg in self.segments:
3258
3259
3260
3261
3262
                try:
                    for d in ('fwd','rev'):
                        seg.strand_pieces[d] = []
                except:
                    pass
3263
        except:
3264
            pass
3265
        self.strands = strands = []
cmaffeo2's avatar
cmaffeo2 committed
3266

3267
3268
3269
        """ Ensure unconnected ends have 5prime Location objects """
        for seg in self.segments:
            ## TODO move into Segment calls
3270
3271
            five_prime_locs = sum([seg.get_locations(s) for s in ("5prime","crossover","terminal_crossover")],[])
            three_prime_locs = sum([seg.get_locations(s) for s in ("3prime","crossover","terminal_crossover")],[])
3272
3273
3274
3275

            def is_start_5prime(l):
                return l.get_nt_pos() < 1 and l.on_fwd_strand
            def is_end_5prime(l):
3276
                return l.get_nt_pos() > seg.num_nt-2 and not l.on_fwd_strand
3277
3278
3279
            def is_start_3prime(l):
                return l.get_nt_pos() < 1 and not l.on_fwd_strand
            def is_end_3prime(l):
3280
                return l.get_nt_pos() > seg.num_nt-2 and l.on_fwd_strand
3281
3282
3283
3284
3285
3286
3287
3288
3289
3290
3291
3292
3293
3294
3295
3296

            if seg.start5.connection is None:
                if len(list(filter( is_start_5prime, five_prime_locs ))) == 0:
                    seg.add_5prime(0) # TODO ensure this is the same place

            if 'end5' in seg.__dict__ and seg.end5.connection is None:
                if len(list(filter( is_end_5prime, five_prime_locs ))) == 0:
                    seg.add_5prime(seg.num_nt-1,on_fwd_strand=False)

            if 'start3' in seg.__dict__ and seg.start3.connection is None:
                if len(list(filter( is_start_3prime, three_prime_locs ))) == 0:
                    seg.add_3prime(0,on_fwd_strand=False)

            if seg.end3.connection is None:
                if len(list(filter( is_end_3prime, three_prime_locs ))) == 0:
                    seg.add_3prime(seg.num_nt-1)
3297
3298
3299
3300
3301
3302
3303
3304

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

cmaffeo2's avatar
cmaffeo2 committed
3305
        """ Build strands from connectivity of helices """
3306
        def _recursively_build_strand(strand, segment, pos, is_fwd, mycounter=0, move_at_least=0.5):
3307
3308
3309
3310
3311
            seg = segment
            history = []
            while True:
                mycounter+=1
                if mycounter > 10000:
3312
                    raise Exception("Too many iterations")
3313
3314
3315
3316
3317
3318
3319

                #if seg.name == "22-1" and pos > 140:
                # if seg.name == "22-2":
                #     import pdb
                #     pdb.set_trace()

                end_pos, next_seg, next_pos, next_dir, move_at_least = seg.get_strand_segment(pos, is_fwd, move_at_least)
3320
                history.append((seg,pos,end_pos,is_fwd))
cmaffeo2's avatar
cmaffeo2 committed
3321
3322
3323
3324
3325
3326
3327
                try:
                    strand.add_dna(seg, pos, end_pos, is_fwd)
                except CircularDnaError:
                    ## Circular DNA was found
                    break
                except:
                    print("Unexpected error:", sys.exc_info()[0])
3328
3329
3330
3331
                    # import pdb
                    # pdb.set_trace()
                    # seg.get_strand_segment(pos, is_fwd, move_at_least)
                    # strand.add_dna(seg, pos, end_pos, is_fwd)
cmaffeo2's avatar
cmaffeo2 committed
3332
                    raise
3333
3334

                if next_seg is None:
3335
                    break
3336
3337
                else:
                    seg,pos,is_fwd = (next_seg, next_pos, next_dir)
3338
            strand.history = list(history)
3339
            return history
3340

3341
3342
        strand_counter = 0
        history = []
cmaffeo2's avatar
cmaffeo2 committed
3343
        for seg in self.segments:
3344
            locs = filter(lambda l: l.connection is None, seg.get_5prime_locations())
3345
3346
3347
            if locs is None: continue
            # for pos, is_fwd in locs:
            for l in locs:
3348
3349
                if _DEBUG_TRACE:
                    print("Tracing 5prime",l)
cmaffeo2's avatar
cmaffeo2 committed
3350
                # TODOTODO
3351
                pos = seg.contour_to_nt_pos(l.address, round_nt=True)
3352
                is_fwd = l.on_fwd_strand
3353
3354
3355
                s = Strand(segname="S{:03d}".format(len(strands)))
                strand_history = _recursively_build_strand(s, seg, pos, is_fwd)
                history.append((l,strand_history))
3356
                # print("{} {}".format(seg.name,s.num_nt))
cmaffeo2's avatar
cmaffeo2 committed
3357
                strands.append(s)
cmaffeo2's avatar
cmaffeo2 committed
3358
3359
3360
3361
3362
3363
3364
3365
3366
3367
3368
3369
3370
3371

        ## Trace circular DNA
        def strands_cover_segment(segment, is_fwd=True):
            direction = 'fwd' if is_fwd else 'rev'
            nt = 0
            for sp in segment.strand_pieces[direction]:
                nt += sp.num_nt
            return nt == segment.num_nt

        def find_nt_not_in_strand(segment, is_fwd=True):
            fwd_str = 'fwd' if is_fwd else 'rev'

            def check(val):
                assert(val >= 0 and val < segment.num_nt)
cmaffeo2's avatar
cmaffeo2 committed
3372
3373
                # print("find_nt_not_in_strand({},{}) returning {}".format(
                #     segment, is_fwd, val))
cmaffeo2's avatar
cmaffeo2 committed
3374
3375
3376
3377
3378
3379
3380
3381
3382
3383
3384
3385
3386
3387
3388
3389
3390
3391
                return val

            if is_fwd:
                last = -1
                for sp in segment.strand_pieces[fwd_str]:
                    if sp.start-last > 1:
                        return check(last+1)
                    last = sp.end
                return check(last+1)
            else:
                last = segment.num_nt
                for sp in segment.strand_pieces[fwd_str]:
                    if last-sp.end > 1:
                        return check(last-1)
                    last = sp.start
                return check(last-1)

        def add_strand_if_needed(seg,is_fwd):
3392
            history = []
cmaffeo2's avatar
cmaffeo2 committed
3393
3394
            if not strands_cover_segment(seg, is_fwd):
                pos = nt = find_nt_not_in_strand(seg, is_fwd)
3395
3396
                if _DEBUG_TRACE:
                    print("Tracing circular strand", pos, is_fwd)
3397
3398
                s = Strand(segname="S{:03d}".format(len(strands)),
                           is_circular = True)
3399
                history = _recursively_build_strand(s, seg, pos, is_fwd)
cmaffeo2's avatar
cmaffeo2 committed
3400
                strands.append(s)
3401
            return history
cmaffeo2's avatar
cmaffeo2 committed
3402
3403
3404
3405
3406
3407
3408
3409
3410
3411
3412
3413
3414
3415
3416
3417
3418
3419
3420

        for seg in self.segments:
            add_strand_if_needed(seg,True)
            if isinstance(seg, DoubleStrandedSegment):
                add_strand_if_needed(seg,False)

        self.strands = sorted(strands, key=lambda s:s.num_nt)[::-1]
        def check_strands():
            dsdna = filter(lambda s: isinstance(s,DoubleStrandedSegment), self.segments)
            for s in dsdna:
                nt_fwd = nt_rev = 0
                for sp in s.strand_pieces['fwd']:
                    nt_fwd += sp.num_nt
                for sp in s.strand_pieces['rev']:
                    nt_rev += sp.num_nt
                assert( nt_fwd == s.num_nt and nt_rev == s.num_nt )
                # print("{}: {},{} (fwd,rev)".format(s.name, nt_fwd/s.num_nt,nt_rev/s.num_nt))
        check_strands()

3421
3422
3423
3424
3425
3426
3427
        ## relabel segname
        counter = 0
        for s in self.strands:
            if s.segname is None:
                s.segname = "D%03d" % counter
            counter += 1

3428
3429
3430
3431
3432
3433
3434
3435
3436
3437
3438
3439
3440
3441
3442
    def _assign_basepairs(self):
        ## Assign basepairs
        for seg in self.segments:
            if isinstance(seg, DoubleStrandedSegment):
                strands1 = seg.strand_pieces['fwd'] # already sorted
                strands2 = seg.strand_pieces['rev']

                nts1 = [nt for s in strands1 for nt in s.children]
                nts2 = [nt for s in strands2 for nt in s.children[::-1]]
                assert(len(nts1) == len(nts2))
                for nt1,nt2 in zip(nts1,nts2):
                    ## TODO weakref
                    nt1.basepair = nt2
                    nt2.basepair = nt1

3443
3444
3445
3446
3447
3448
3449
3450
3451
3452
3453
3454
3455
3456
3457
3458
3459
3460
3461
3462
3463
3464
3465
3466
3467
3468
3469
    def write_atomic_bp_restraints(self, output_name, spring_constant=1.0):
        ## TODO: ensure atomic model was generated already
        ## TODO: allow ENM to be created without first building atomic model

        with open("%s.exb" % output_name,'w') as fh:
            for seg in self.segments:
                ## Continue unless dsDNA
                if not isinstance(seg,DoubleStrandedSegment): continue
                for strand_piece in seg.strand_pieces['fwd']:
                    assert( strand_piece.is_fwd )
                    for nt1 in strand_piece.children:
                        nt2 = nt1.basepair
                        if nt1.resname == 'ADE':
                            names = (('N1','N3'),('N6','O4'))
                        elif nt1.resname == 'GUA':
                            names = (('N2','O2'),('N1','N3'),('O6','N4'))
                        elif nt1.resname == 'CYT':
                            names = (('O2','N2'),('N3','N1'),('N4','O6'))
                        elif nt1.resname == 'THY':
                            names = (('N3','N1'),('O4','N6'))
                        else:
                            raise Exception("Unrecognized nucleotide!")
                        for n1, n2 in names:
                            i = nt1._get_atomic_index(name=n1)
                            j = nt2._get_atomic_index(name=n2)
                            fh.write("bond %d %d %f %.2f\n" % (i,j,spring_constant,2.8))

3470
    def write_atomic_ENM(self, output_name, lattice_type=None, interhelical_bonds=True):
3471
3472
        ## TODO: ensure atomic model was generated already
        if lattice_type is None:
3473
3474
3475
3476
            try:
                lattice_type = self.lattice_type
            except:
                lattice_type = "square"
cmaffeo2's avatar
cmaffeo2 committed
3477
3478
3479
3480
3481
3482
        else:
            try:
                if lattice_type != self.lattice_type:
                    print("WARNING: printing ENM with a lattice type ({}) that differs from model's lattice type ({})".format(lattice_type,self.lattice_type))
            except:
                pass
3483
3484
3485
3486
3487
3488
3489
3490

        if lattice_type == "square":
            enmTemplate = enmTemplateSQ
        elif lattice_type == "honeycomb":
            enmTemplate = enmTemplateHC
        else:
            raise Exception("Lattice type '%s' not supported" % self.latticeType)

3491
        ## TODO: allow ENM to be created without first building atomic model
3492
3493
        noStackPrime = 0
        noBasepair = 0
3494
        with open("%s.exb" % output_name,'w') as fh:
3495
3496
            # natoms=0

3497
3498
3499
3500
3501
3502
3503
3504
            for seg in self.segments:
                ## Continue unless dsDNA
                if not isinstance(seg,DoubleStrandedSegment): continue
                for strand_piece in seg.strand_pieces['fwd'] + seg.strand_pieces['rev']:
                    for nt1 in strand_piece.children:
                        other = []
                        nt2 = nt1.basepair
                        if strand_piece.is_fwd:
3505
                            other.append((nt2,'pair'))
3506
3507
3508
3509

                        nt2 = nt2.get_intrahelical_above()
                        if nt2 is not None and strand_piece.is_fwd:
                            ## TODO: check if this already exists
3510
3511
                            other.append((nt2,'paircross'))

3512
3513
3514
                        nt2 = nt1.get_intrahelical_above()
                        if nt2 is not None:
                            other.append((nt2,'stack'))
3515
3516
3517
3518
3519
3520
                            try:
                                nt2 = nt2.basepair
                            except:
                                assert( isinstance(nt2.parent.segment, SingleStrandedSegment) )
                                nt2 = None

3521
                            if nt2 is not None and strand_piece.is_fwd:
3522
3523
3524
3525
3526
3527
3528
3529
3530
3531
3532
3533
3534
3535
3536
                                other.append((nt2,'cross'))

                        for nt2,key in other:
                            """
                            if np.linalg.norm(nt2.position-nt1.position) > 7:
                                import pdb
                                pdb.set_trace()
                            """
                            key = ','.join((key,nt1.sequence[0],nt2.sequence[0]))
                            for n1, n2, d in enmTemplate[key]:
                                d = float(d)
                                k = 0.1
                                if lattice_type == 'honeycomb':
                                    correctionKey = ','.join((key,n1,n2))
                                    assert(correctionKey in enmCorrectionsHC)
3537
3538
3539
                                    dk,dr = enmCorrectionsHC[correctionKey]
                                    k  = float(dk)
                                    d += float(dr)
3540

3541
3542
3543
3544
3545
3546
3547
3548
3549
3550
3551
3552
3553
3554
3555
3556
3557
3558
3559
3560
3561
3562
3563
3564
3565
3566
3567
3568
3569
3570
3571
3572
3573
3574
3575
3576
                                i = nt1._get_atomic_index(name=n1)
                                j = nt2._get_atomic_index(name=n2)
                                fh.write("bond %d %d %f %.2f\n" % (i,j,k,d))

            # print("NO STACKS found for:", noStackPrime)
            # print("NO BASEPAIRS found for:", noBasepair)

        ## Loop dsDNA regions
        push_bonds = []
        processed_segs = set()
        ## TODO possibly merge some of this code with SegmentModel.get_consecutive_crossovers()
        for segI in self.segments: # TODOTODO: generalize through some abstract intrahelical interface that effectively joins "segments", for now interhelical bonds that cross intrahelically-connected segments are ignored
            if not isinstance(segI,DoubleStrandedSegment): continue

            ## Loop over dsDNA regions connected by crossovers
            conn_locs = segI.get_contour_sorted_connections_and_locations("crossover")
            other_segs = list(set([B.container for c,A,B in conn_locs]))
            for segJ in other_segs:
                if (segI,segJ) in processed_segs:
                    continue
                processed_segs.add((segI,segJ))
                processed_segs.add((segJ,segI))

                ## TODO perhaps handle ends that are not between crossovers

                ## Loop over ordered pairs of crossovers between the two
                cls = filter(lambda x: x[-1].container == segJ, conn_locs)
                cls = sorted( cls, key=lambda x: x[1].get_nt_pos() )
                for cl1,cl2 in zip(cls[:-1],cls[1:]):
                    c1,A1,B1 = cl1
                    c2,A2,B2 = cl2

                    ntsI1,ntsI2 = [segI.contour_to_nt_pos(A.address) for A in (A1,A2)]
                    ntsJ1,ntsJ2 = [segJ.contour_to_nt_pos(B.address) for B in (B1,B2)]
                    ntsI = ntsI2-ntsI1+1
                    ntsJ = ntsJ2-ntsJ1+1
3577
3578
                    # assert( np.isclose( ntsI, int(round(ntsI)) ) )
                    # assert( np.isclose( ntsJ, int(round(ntsJ)) ) )
3579
3580
3581
3582
3583
3584
3585
3586
3587
3588
3589
3590
3591
3592
3593
3594
                    ntsI,ntsJ = [int(round(i)) for i in (ntsI,ntsJ)]

                    ## Find if dsDNA "segments" are pointing in same direction
                    ## could move this block out of the loop
                    tangentA = segI.contour_to_tangent(A1.address)
                    tangentB = segJ.contour_to_tangent(B1.address)
                    dot1 = tangentA.dot(tangentB)
                    tangentA = segI.contour_to_tangent(A2.address)
                    tangentB = segJ.contour_to_tangent(B2.address)
                    dot2 = tangentA.dot(tangentB)

                    if dot1 > 0.5 and dot2 > 0.5:
                        ...
                    elif dot1 < -0.5 and dot2 < -0.5:
                        ## TODO, reverse
                        ...
3595
                        # print("Warning: {} and {} are on antiparallel helices (not yet implemented)... skipping".format(A1,B1))
3596
3597
                        continue
                    else:
3598
                        # print("Warning: {} and {} are on helices that do not point in similar direction... skipping".format(A1,B1))
3599
3600
                        continue

3601
3602
3603
3604
3605
3606
3607
3608
3609
3610
3611
3612
3613
3614
3615
3616
3617
3618
3619
3620
3621
3622
3623
3624
                    ## Go through each nucleotide between the two
                    for ijmin in range(min(ntsI,ntsJ)):
                        i=j=ijmin
                        if ntsI < ntsJ:
                            j = int(round(float(ntsJ*i)/ntsI))
                        elif ntsJ < ntsI:
                            i = int(round(float(ntsI*j)/ntsJ))
                        ntI_idx = int(round(ntsI1+i))
                        ntJ_idx = int(round(ntsJ1+j))

                        ## Skip nucleotides that are too close to crossovers
                        if i < 11 or j < 11: continue
                        if ntsI2-ntI_idx < 11 or ntsJ2-ntJ_idx < 11: continue

                        ## Find phosphates at ntI/ntJ
                        for direction in [True,False]:
                            try:
                                i = segI._get_atomic_nucleotide(ntI_idx, direction)._get_atomic_index(name="P")
                                j = segJ._get_atomic_nucleotide(ntJ_idx, direction)._get_atomic_index(name="P")
                                push_bonds.append((i,j))
                            except:
                                # print("WARNING: could not find 'P' atom in {}:{} or {}:{}".format( segI, ntI_idx, segJ, ntJ_idx ))
                                ...

3625
        # print("PUSH BONDS:", len(push_bonds))
3626
3627
3628
3629
3630
3631
3632
3633
3634
3635
3636
3637
3638
3639
        if interhelical_bonds:
            if not self.useTclForces:
                with open("%s.exb" % output_name, 'a') as fh:
                    for i,j in push_bonds:
                        fh.write("bond %d %d %f %.2f\n" % (i,j,1.0,31))
            else:
                flat_push_bonds = list(sum(push_bonds))
                atomList = list(set( flat_push_bonds ))
                with open("%s.forces.tcl" % output_name,'w') as fh:
                    fh.write("set atomList {%s}\n\n" %
                             " ".join([str(x-1) for x in  atomList]) )
                    fh.write("set bonds {%s}\n" %
                             " ".join([str(x-1) for x in flat_push_bonds]) )
                    fh.write("""
3640
3641
3642
3643
3644
3645
3646
3647
3648
3649
3650
3651
3652
3653
3654
3655
3656
3657
3658
3659
3660
3661
3662
3663
3664
3665
3666
3667
3668
foreach atom $atomList {
    addatom $atom
}

proc calcforces {} {
    global atomList bonds
    loadcoords rv

    foreach i $atomList {
        set force($i) {0 0 0}
    }

    foreach {i j} $bonds {
        set rvec [vecsub $rv($j) $rv($i)]
        # lassign $rvec x y z
        # set r [expr {sqrt($x*$x+$y*$y+$z*$z)}]
        set r [getbond $rv($j) $rv($i)]
        set f [expr {2*($r-31.0)/$r}]
        vecadd $force($i) [vecscale $f $rvec]
        vecadd $force($j) [vecscale [expr {-1.0*$f}] $rvec]
    }

    foreach i $atomList {
        addforce $i $force($i)
    }

}
""")

3669
3670
3671
    def get_bounding_box( self, num_points=3 ):
        positions = np.zeros( (len(self.segments)*num_points, 3) )
        i = 0
3672
        for s in self.segments:
3673
3674
3675
3676
3677
3678
3679
3680
3681
3682
3683
3684
3685
3686
            for c in np.linspace(0,1,num_points):
                positions[i] = (s.contour_to_position(c))
                i += 1
        min_ = np.array([np.min(positions[:,i]) for i in range(3)])
        max_ = np.array([np.max(positions[:,i]) for i in range(3)])
        return min_,max_

    def get_bounding_box_center( self, num_points=3 ):
        min_,max_ = self.get_bounding_box(num_points)
        return 0.5*(max_+min_)

    def dimensions_from_structure( self, padding_factor=1.5, isotropic=False ):
        min_,max_ = self.get_bounding_box()
        dx,dy,dz = (max_-min_+30)*padding_factor
3687
3688
        if isotropic:
            dx = dy = dz = max((dx,dy,dz))
3689
        return np.array([dx,dy,dz])
3690

3691
    def add_grid_potential(self, grid_file, scale=1, per_nucleotide=True, filter_fn=None):
3692
3693
3694
3695
3696
        grid_file = Path(grid_file)
        if not grid_file.is_file():
            raise ValueError("Grid file {} does not exist".format(grid_file))
        if not grid_file.is_absolute():
            grid_file = Path.cwd() / grid_file
3697
        self.grid_potentials.append((grid_file,scale,per_nucleotide, filter_fn))
3698

3699
    def _apply_grid_potentials_to_beads(self, bead_type_dict ):
3700
3701
3702
        if len(self.grid_potentials) > 1:
            raise NotImplementedError("Multiple grid potentials are not yet supported")

3703
3704
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
3721
3722
3723
3724
3725
        def add_grid_to_type(particle_type):
            if particle_type.name[0] == "O": return
            s = scale*particle_type.nts if per_nucleotide else scale
            try:
                particle_type.grid = particle_type.grid + (grid_file, s)
            except:
                particle_type.grid = tuple((grid_file, s))

        for grid_file, scale, per_nucleotide, filter_fn in self.grid_potentials:
            if filter_fn is None:
                for key,particle_type in bead_type_dict.items():
                    add_grid_to_type(particle_type)
            else:
                grid_types = dict()
                for b in filter(filter_fn, sum([seg.beads for seg in self.segments],[])):
                    t = b.type_
                    if t.name[0] == "O": continue
                    if t not in grid_types:
                        new_type = ParticleType(name=t.name+'G',charge=t.charge, parent=t)
                        add_grid_to_type(new_type)
                        grid_types[t] = new_type
                    b.type_ = grid_types[t]

3726

3727
    def _generate_atomic_model(self, scale=1):
3728
3729
3730
3731
3732
3733
        ## TODO: deprecate
        self.generate_atomic_model(scale=scale)

    def generate_atomic_model(self, scale=1):
        self.clear_beads()
        self.children = self.strands # TODO: is this going to be okay? probably
3734
        first_atomic_index = 0
3735
        for s in self.strands:
3736
3737
            first_atomic_index = s.generate_atomic_model(scale,first_atomic_index)
        self._assign_basepairs()
3738

3739
3740
    """ OxDNA """
    ## https://dna.physics.ox.ac.uk/index.php/Documentation#Input_file
3741
3742
3743
3744
3745
3746
    def generate_oxdna_model(self, scale=1):
        self.clear_beads()
        self.children = self.strands
        for s in self.strands:
            s.generate_oxdna_model()

3747
3748
3749
3750
3751
3752
3753
3754
3755
3756
3757
3758
3759
3760
3761
3762
3763
3764
3765
3766
3767
3768
3769
3770
3771
3772
3773
3774
3775
3776
3777
3778
3779
3780
3781
3782
3783
3784
3785
3786
3787
    def _write_oxdna_configuration(self, filename):

        _angstroms_to_oxdna = 0.11739845 ##  units "AA" "8.518e-10 m"
        with open(filename,'w') as fh:
            fh.write("""t = {temperature}
b = {dimX} {dimY} {dimZ}
E = 0 0 0
""".format(temperature=self.temperature,
                        dimX = self.dimensions[0]*_angstroms_to_oxdna,
                        dimY = self.dimensions[1]*_angstroms_to_oxdna,
                        dimZ = self.dimensions[2]*_angstroms_to_oxdna))

            base_vec = np.array((1,0,0)) # TODO
            norm_vec = np.array((0,0,1)) # TODO

            ## Traverse 3'-to-5'
            for nt in [nt for s in self.strands for nt in s.oxdna_nt[::-1]]:
                data = dict()
                # o = nt.collapsedOrientation()
                o = nt.orientation
                for k,v in zip('x y z'.split(),nt.collapsedPosition()):
                    data[k] = v * _angstroms_to_oxdna
                for k,v in zip('x y z'.split(),o.dot(base_vec)):
                    data['b'+k] = v
                for k,v in zip('x y z'.split(),o.dot(norm_vec)):
                    data['n'+k] = v
                fh.write("{x} {y} {z} {bx} {by} {bz} {nx} {ny} {nz}   0 0 0   0 0 0\n".format(**data)
    )

    def _write_oxdna_topology(self,filename):

        strands = [s for s in self.strands if s.num_nt > 0]
        with open(filename,'w') as fh:

            fh.write("{num_nts} {num_strands}\n".format(
                num_nts = sum([s.num_nt for s in strands]),
                num_strands = len(self.strands)))

            idx = 0
            sidx = 1
            for strand in strands:
3788
                prev = idx+strand.num_nt-1 if strand.is_circular else -1
3789
3790
3791
3792
3793
3794
3795
3796
3797
3798
3799
3800
3801
3802
3803
3804
3805
3806
3807
3808
3809
3810
3811
3812
3813
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
3826
3827
                last = idx if strand.is_circular else -1

                ## Traverse 3'-to-5'
                sequence = [seq for s in strand.strand_segments
                            for seq in s.get_sequence()][::-1]
                for seq in sequence[:-1]:
                    ## strand seq 3' 5'
                    fh.write("{} {} {} {}\n".format(sidx, seq, prev, idx+1))
                    prev = idx
                    idx += 1
                seq = sequence[-1]
                fh.write("{} {} {} {}\n".format(sidx, seq, prev, last))
                idx += 1
                sidx += 1

    def _write_oxdna_input(self, filename,
                           topology,
                           conf_file,
                           trajectory_file,
                           last_conf_file,
                           log_file,
                           num_steps = 1e6,
                           interaction_type = 'DNA2',
                           salt_concentration = None,
                           print_conf_interval = None,
                           print_energy_every = None,
                           timestep = 0.003,
                           sim_type = "MD",
                           backend = None,
                           backend_precision = None,
                           seed = None,
                           newtonian_steps = 103,
                           diff_coeff = 2.50,
                           thermostat = "john",
                           list_type = "cells",
                           ensemble = "nvt",
                           delta_translation = 0.22,
                           delta_rotation = 0.22,
                           verlet_skin = 0.5,
3828
                           max_backbone_force = 100,
3829
3830
                           external_forces_file = None,
                           seq_dep_file = None
3831
3832
3833
3834
3835
3836
3837
3838
3839
3840
3841
3842
3843
3844
3845
3846
3847
3848
3849
3850
3851
3852
3853
3854
3855
3856
3857
3858
3859
3860
3861
3862
3863
3864
3865
3866
3867
3868
3869
3870
3871
3872
3873
3874
3875
3876
3877
3878
3879
                       ):

        if seed is None:
            import random
            seed = random.randint(1,99999)

        temperature = self.temperature
        num_steps = int(num_steps)
        newtonian_steps = int(newtonian_steps)

        if print_conf_interval is None:
            # print_conf_interval = min(num_steps//100)
            print_conf_interval = 10000
        print_conf_interval = int(print_conf_interval)

        if print_energy_every is None:
            print_energy_every = print_conf_interval
        print_energy_every = int(print_energy_every)

        if max_backbone_force is not None:
            max_backbone_force = 'max_backbone_force = {}'.format(max_backbone_force)


        if interaction_type in ('DNA2',):
            if salt_concentration is None:
                try:
                    ## units "80 epsilon0 295 k K / (2 (AA)**2 e**2/particle)" mM
                    salt_concentration = 9331.3126/self.debye_length**2
                except:
                    salt_concentration = 0.5
            salt_concentration = 'salt_concentration = {}'.format(salt_concentration)
        else:
            salt_concentration = ''

        if backend is None:
            backend = 'CUDA' if sim_type == 'MD' else 'CPU'

        if backend_precision is None:
            backend_precision = 'mixed' if backend == 'CUDA' else 'double'

        if sim_type == 'VMMC':
            ensemble = 'ensemble = {}'.format(ensemble)
            delta_translation = 'delta_translation = {}'.format(delta_translation)
            delta_rotation = 'delta_rotation = {}'.format(delta_rotation)
        else:
            ensemble = ''
            delta_translation = ''
            delta_rotation = ''

3880
3881
3882
3883
3884
        if external_forces_file is None:
            external_forces = "external_forces = 0"
        else:
            external_forces = "external_forces = 1\nexternal_forces_file = {}".format(external_forces_file)

3885
3886
3887
3888
3889
3890
3891
3892
3893
3894
        if seq_dep_file is None:
            sequence_dependence = ""
        else:

            if seq_dep_file == "oxDNA2":
                seq_dep_file = get_resource_path("oxDNA2_sequence_dependent_parameters.txt")
            elif seq_dep_file == "oxDNA1":
                seq_dep_file = get_resource_path("oxDNA1_sequence_dependent_parameters.txt")
            sequence_dependence = "use_average_seq = 1\nseq_dep_file = {}".format(seq_dep_file)

3895
3896
3897
3898
3899
3900
3901
3902
3903
3904
3905
3906
3907
3908
3909
3910
3911
3912
3913
3914
3915
3916
3917
3918
3919
3920
3921
3922
3923
3924
        with open(filename,'w') as fh:
            fh.write("""##############################
####  PROGRAM PARAMETERS  ####
##############################
interaction_type = {interaction_type}
{salt_concentration}
sim_type = {sim_type}
backend = {backend}
backend_precision = {backend_precision}
#debug = 1
seed = {seed}

##############################
####    SIM PARAMETERS    ####
##############################
steps = {num_steps:d}
newtonian_steps = {newtonian_steps:d}
diff_coeff = {diff_coeff}
thermostat = {thermostat}

list_type = {list_type}
{ensemble}
{delta_translation}
{delta_rotation}

T = {temperature:f} K
dt = {timestep}
verlet_skin = {verlet_skin}
{max_backbone_force}

3925
{external_forces}
3926
{sequence_dependence}
3927

3928
3929
3930
3931
3932
3933
3934
3935
3936
3937
3938
3939
3940
3941
3942
3943
3944
3945
3946
3947
3948
3949
3950
3951
3952
3953
3954
3955
3956
3957
3958
3959
3960
3961
3962
3963
3964
3965
3966
3967
3968
3969
3970
3971
3972
3973
3974
3975
3976
3977
3978
3979
3980
3981
3982
3983
3984
3985
3986
3987
3988
3989
3990
3991
3992
3993
3994
3995
3996
3997
3998
3999
##############################
####    INPUT / OUTPUT    ####
##############################
topology = {topology}
conf_file = {conf_file}
lastconf_file = {last_conf_file}
trajectory_file = {trajectory_file}
refresh_vel = 1
log_file = {log_file}
no_stdout_energy = 1
restart_step_counter = 1
energy_file = {log_file}.energy.dat
print_conf_interval = {print_conf_interval}
print_energy_every = {print_energy_every}
time_scale = linear
""".format( **locals() ))

    def simulate_oxdna(self, output_name, directory='.', output_directory='output', topology=None, configuration=None, oxDNA=None, **oxdna_args):

        if output_directory == '': output_directory='.'
        d_orig = os.getcwd()
        if not os.path.exists(directory):
            os.makedirs(directory)

        os.chdir(directory)
        try:

            if oxDNA is None:
                for path in os.environ["PATH"].split(os.pathsep):
                    path = path.strip('"')
                    fname = os.path.join(path, "oxDNA")
                    if os.path.isfile(fname) and os.access(fname, os.X_OK):
                        oxDNA = fname
                        break

            if oxDNA is None: raise Exception("oxDNA was not found")

            if not os.path.exists(oxDNA):
                raise Exception("oxDNA was not found")
            if not os.path.isfile(oxDNA):
                raise Exception("oxDNA was not found")
            if not os.access(oxDNA, os.X_OK):
                raise Exception("oxDNA is not executable")

            if not os.path.exists(output_directory):
                os.makedirs(output_directory)
            elif not os.path.isdir(output_directory):
                raise Exception("output_directory '%s' is not a directory!" % output_directory)


            if configuration is None:
                configuration = "{}.conf".format(output_name)
                self._write_oxdna_configuration(configuration)
            # elif not Path(configuration).exists():
            #     raise Exception("Unable to find oxDNA configuration file '{}'.".format(configuration))

            if topology is None:
                topology = "{}-topology.dat".format(output_name)
                self._write_oxdna_topology(topology)
            elif not Path(topology).exists():
                raise Exception("Unable to find oxDNA topology file '{}'.".format(topology))

            last_conf_file = "{}/{}.last.conf".format(output_directory,output_name)
            input_file = "{}-input".format(output_name)

            self._write_oxdna_input(input_file,
                                    topology = topology,
                                    conf_file = configuration,
                                    trajectory_file = "{}/{}.dat".format(output_directory,output_name),
                                    last_conf_file = last_conf_file,
                                    log_file="{}/{}.log".format(output_directory,output_name),
                                    **oxdna_args)
4000
            # os.sync()
For faster browsing, not all history is shown. View entire blame