Fix multiple issues related to strand generation, also reduce timestep to...

Fix multiple issues related to strand generation, also reduce timestep to 100fs if coarse_local_twist==True

mrdna/arbdmodel/__init__.py

@@ -8,6 +8,7 @@ import numpy as np
 from copy import copy, deepcopy
 from inspect import ismethod
 import os, sys, subprocess
+import shutil
 
 _RESOURCE_DIR = Path(__file__).parent / 'resources'
 def get_resource_path(relative_path):
@@ -416,6 +417,28 @@ class ParticleType():
         else:
             return False
 
+    def __getattr__(self, name):
+        """
+        Try to get attribute from the parent
+
+        """
+        if "parent" not in self.__dict__ or self.__dict__["parent"] is None or name == "children":
+            raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name))
+
+        excluded_attributes = ParticleType.excludedAttributes
+        if name in excluded_attributes:
+            raise AttributeError("'{}' object has no attribute '{}' and cannot look it up from the parent".format(type(self).__name__, name))
+
+        ## TODO: determine if there is a way to avoid __getattr__ if a method is being looked up
+        try:
+            ret = getattr(self.parent,name)
+        except:
+            raise AttributeError("'{}' object has no attribute '{}'".format(type(self).__name__, name))
+        if ismethod(ret):
+            raise AttributeError("'{}' object has no method '{}'".format(type(self).__name__, name))
+        return ret 
+
+        
     def __hash_key(self):
         l = [self.name,self.charge]
         for keyval in sorted(self.__dict__.items()):
@@ -1463,7 +1486,7 @@ if {{$nLast == 0}} {{
 run {num_steps:d}
 """.format(**format_data))
 
-    def atomic_simulate(self, output_name, output_directory='output', dry_run = False, namd2=None, log_file=None, num_procs=None, gpu=None):
+    def atomic_simulate(self, output_name, output_directory='output', dry_run = False, namd2=None, log_file=None, num_procs=None, gpu=None, copy_ff_from=get_resource_path("charmm36.nbfix") ):
         if self.cacheUpToDate == False: # TODO: remove cache?
             self._countParticleTypes()
             self._updateParticleOrder()
@@ -1473,6 +1496,14 @@ run {num_steps:d}
         self.writePdb( output_name + ".fixed.pdb", beta_from_fixed=True )
         self.writePsf( output_name + ".psf" )
         self.write_namd_configuration( output_name, output_directory = output_directory )
+
+        if copy_ff_from is not None and copy_ff_from != '':
+            try:
+                shutil.copytree( copy_ff_from, Path(copy_ff_from).stem )
+            except FileExistsError:
+                pass
+
+        
         # os.sync()
 
         if not dry_run:

mrdna/readers/segmentmodel_from_pdb.py

@@ -112,6 +112,10 @@ def find_base_position_orientation( u, selection_text='all' ):
         assert( (ref.atoms.names == sel.atoms.names[sel_order]).all() )
 
         ## Reorder coordinates in sel so atom names match 
+        if np.all(sel.atoms.masses == 0):
+            sel.atoms.masses = [1]*len(sel.atoms)
+        if np.all(ref.atoms.masses == 0):
+            ref.atoms.masses = [1]*len(ref.atoms)
         c = sel.atoms.center_of_mass()
         R,rmsd = align.rotation_matrix(sel.positions[sel_order] - c,
                                        ref.positions - ref.atoms.center_of_mass())                                

mrdna/segmentmodel.py

@@ -100,6 +100,7 @@ class Location():
     def delete(self):
         if self.connection is not None:
             self.connection.delete()
+            self.connection = None
         if self.container is not None:
             self.container.locations.remove(self)
 
@@ -136,7 +137,10 @@ class Connection():
             self.B.container.connections.remove(self)
         self.A.connection = None
         self.B.connection = None
-
+        self.A.type_ = "3prime" # TODO move specialization to inherited class
+        self.B.type_ = "5prime"
+        self.A = self.B = None
+        
     def __repr__(self):
         return "<Connection {}--{}--{}]>".format( self.A, self.type_, self.B )
         
@@ -163,7 +167,7 @@ class ConnectableElement():
                 counter[l] = 1
         assert( np.all( [counter[l] == 1 for l in locs] ) )
         return locs
-
+        
     def get_location_at(self, address, on_fwd_strand=True, new_type=None):
         loc = None
         if (self.num_nt == 1):
@@ -502,6 +506,10 @@ class Segment(ConnectableElement, Group):
         max_ = np.array([np.max(positions[:,i]) for i in range(3)])
         return min_,max_
 
+    def get_location_at_nt(self, nt, on_fwd_strand=True, new_type=None):
+        a = self.nt_pos_to_contour(nt)
+        return self.get_location_at(a, on_fwd_strand, new_type)
+        
     def _get_location_positions(self):
         return [self.contour_to_nt_pos(l.address) for l in self.locations]
 
@@ -547,13 +555,18 @@ class Segment(ConnectableElement, Group):
         if first_nt == last_nt:
             return
 
-        def _transform_contour_positions(contours):
+        def _transform_contour_positions(contours, locations = None):
             nt = self.contour_to_nt_pos(contours)
             first = first_nt if first_nt > 0 else -np.inf
             last = last_nt if last_nt < self.num_nt-1 else np.inf
             bad_ids = (nt >= first-0.5) * (nt <= last+0.5)
             nt[nt>last] = nt[nt>last]-removed_nt
             nt[bad_ids] = np.nan
+            # if locations is not None:
+            #     for i in np.where(nt == 0):
+            #         l = locations[i]
+            #         l.
+                    
             return self.nt_pos_to_contour(nt)* self.num_nt/num_nt
 
         first_nt = np.around(first_nt)
@@ -711,8 +724,14 @@ class Segment(ConnectableElement, Group):
     def contour_to_nt_pos(self, contour_pos, round_nt=False):
         nt = contour_pos*(self.num_nt) - 0.5
         if round_nt:
-            assert( np.isclose(np.around(nt),nt) and contour_pos not in (0,1)  )
-            nt = np.around(nt)
+            # assert( np.isclose(np.around(nt),nt) and contour_pos not in (0,1)  )
+            if contour_pos == 0:
+                nt = 0
+            elif contour_pos == 1:
+                nt = self.num_nt-1
+            else:
+                assert( np.isclose(np.around(nt),nt) )
+                nt = np.around(nt)
         return nt
 
     def nt_pos_to_contour(self,nt_pos):
@@ -944,7 +963,7 @@ class Segment(ConnectableElement, Group):
         ## validate the input
         for end in (end1, end2):
             assert( isinstance(end, Location) )
-            assert( end.type_ in ("end3","end5") )
+            assert( end.type_ in ("end3","end5",'3prime','5prime') )
         assert( end1.type_ != end2.type_ )
 
         if ((isinstance(end1.container, SingleStrandedSegment) and
@@ -1560,6 +1579,11 @@ class StrandInSegment(Group):
         self.end = end
         self.is_fwd = is_fwd
 
+        if is_fwd:
+            assert(end>=start)
+        else:
+            assert(start>=end)
+        
         nts = np.abs(end-start)+1
         self.num_nt = int(round(nts))
         assert( np.isclose(self.num_nt,nts) )
@@ -1826,6 +1850,9 @@ class SegmentModel(ArbdModel):
 
         # self.max_basepairs_per_bead = max_basepairs_per_bead     # dsDNA
         # self.max_nucleotides_per_bead = max_nucleotides_per_bead # ssDNA
+        if isinstance(segments,Segment):
+            segments = [segments]
+
         self.children = self.segments = segments
 
         if (hj_equilibrium_angle > 180) or (hj_equilibrium_angle < -180):
@@ -1967,7 +1994,10 @@ class SegmentModel(ArbdModel):
 
         return twist_spring_from_lp(sep, Lp, self.temperature)
 
-    def extend(self, other, copy=True, include_strands=True):
+    def extend(self, other, copy=False, include_strands=True):
+        if copy == True:
+            logger.warning("Forcing SegmentModel.extend(...,copy=False,...)")
+            copy = False
         assert( isinstance(other, SegmentModel) )
 
         try:
@@ -2325,39 +2355,73 @@ class SegmentModel(ArbdModel):
 
     def convert_crossover_to_intrahelical(self, crossover):
         c = crossover
-        seg1 = c.A.container
-        seg2 = c.B.container
-        nt1 = c.A.get_nt_pos()
-        nt2 = c.B.get_nt_pos()
-
-        assert( c.A.on_fwd_strand == c.B.on_fwd_strand )
-        on_fwd_strand = c.A.on_fwd_strand
+        A = c.A if c.B.is_3prime_side_of_connection else c.B
+        B = c.B if c.B.is_3prime_side_of_connection else c.A
+        assert( not A.is_3prime_side_of_connection )
+        assert(     B.is_3prime_side_of_connection )
+        
+        seg1 = A.container
+        seg2 = B.container
+        nt1 = A.get_nt_pos()
+        nt2 = B.get_nt_pos()
 
-        c.delete()
+        # assert( c.A.on_fwd_strand == c.B.on_fwd_strand )
+        on_fwd_strand_A = A.on_fwd_strand
+        on_fwd_strand_B = B.on_fwd_strand
 
-        if nt1 < 1 and nt2 > seg2.num_nt-2:
-            if on_fwd_strand:
-                seg2.connect_end3(seg1.start5)
+        # print(nt1, nt2, on_fwd_strand_A, on_fwd_strand_B, c)
+        
+        if on_fwd_strand_A:
+            assert( nt1 > seg1.num_nt-2 )
+            fn = seg1.connect_end3
+            if on_fwd_strand_B:
+                assert( nt2 < 1 )
+                arg = seg2.start5
             else:
-                seg1.connect_start3(seg2.end5)
-
-        elif nt2 < 1 and nt1 > seg1.num_nt-2:
-            if on_fwd_strand:
-                seg1.connect_end3(seg2.start5)
+                assert( nt2 > seg2.num_nt-2 )                
+                arg = seg2.end5
+        else:
+            assert( nt1 < 1 )
+            fn = seg1.connect_start3
+            if on_fwd_strand_B:
+                assert( nt2 < 1 )
+                arg = seg2.start5
             else:
-                seg2.connect_start3(seg1.end5)
+                assert( nt2 > seg2.num_nt-2 )                
+                arg = seg2.end5
+        c.delete()
+        fn(arg)
 
     def get_blunt_DNA_ends(self):
         return sum([s.get_blunt_DNA_ends for s in self.segments],[])
 
     def convert_crossovers_at_ends_beyond_cutoff_to_intrahelical(self,cutoff):
+        if cutoff < 0:
+            raise ValueError("Cutoff should be positive!")
+        return convert_crossovers_to_intrahelical(
+            position_filter = lambda r1,r2: np.linalg.norm(r2-r1) > cutoff,
+            terminal_only = True)
+
+
+    def convert_crossovers_to_intrahelical(self, position_filter=None, terminal_only=False):
+        conn = self.get_crossovers_at_ends() if terminal_only else [c for c,A,B in self.get_connections('crossover')]
+        for c in conn:
+            if position_filter is not None:
+                r1,r2 = [l.container.contour_to_position(l.address)
+                         for l in (c.A,c.B)]
+                if position_filter(r1,r2):
+                    self.convert_crossover_to_intrahelical(c)
+            else:
+                self.convert_crossover_to_intrahelical(c)
+
+    def convert_close_crossovers_to_intrahelical(self,cutoff):
         if cutoff < 0:
             raise ValueError("Cutoff should be positive!")
 
-        for c in self.get_crossovers_at_ends():
+        for c,A,B in self.get_connections('crossover'):
             r1,r2 = [l.container.contour_to_position(l.address)
                      for l in (c.A,c.B)]
-            if np.linalg.norm(r2-r1) > cutoff:
+            if np.linalg.norm(r2-r1) < cutoff:
                 self.convert_crossover_to_intrahelical(c)
 
     def _generate_bead_model(self,
@@ -2556,7 +2620,7 @@ class SegmentModel(ArbdModel):
                 ## already handled by Segment._generate_beads, unless A,B are on same segment
                 if A.container == B.container:
                     sorted_sites = sorted([(abs(L.address-b1.contour_position),L)
-                                       for L in (A,B)])
+                                           for L in (A,B)], key=lambda x: x[0])
                     close_site, far_site = [s[1] for s in sorted_sites]
 
                     b3 = far_site.container.get_nearest_bead(far_site.address)
@@ -2804,6 +2868,11 @@ class SegmentModel(ArbdModel):
         beads = dists.keys()
         def _recursively_get_beads_within(b1,d,done=()):
             ret = []
+            self
+            intra_beads
+            if b1 not in dists:
+                print("Skipping exclusions for",b1)
+                return ret
             for b2,sep in dists[b1].items():
                 if b2 in done: continue
                 if sep < d:
@@ -3184,6 +3253,10 @@ class SegmentModel(ArbdModel):
                     pot = self.get_dihedral_potential(k,t0)
                     self.add_dihedral( n1,n2,n3,n4, pot )
 
+        for seg in self.segments:
+            for callback in seg._generate_bead_callbacks:
+                callback(seg)
+                    
         for callback in self._generate_bead_callbacks:
             callback(self)
 
@@ -3330,13 +3403,17 @@ class SegmentModel(ArbdModel):
             three_prime_locs = sum([seg.get_locations(s) for s in ("3prime","crossover","terminal_crossover")],[])
 
             def is_start_5prime(l):
-                return l.get_nt_pos() < 1 and l.on_fwd_strand
+                is_end =  l.get_nt_pos() < 1 and l.on_fwd_strand
+                return is_end and (l.connection is None or l.is_3prime_side_of_connection is not False)
             def is_end_5prime(l):
-                return l.get_nt_pos() > seg.num_nt-2 and not l.on_fwd_strand
+                is_end = l.get_nt_pos() > seg.num_nt-2 and not l.on_fwd_strand
+                return is_end and (l.connection is None or l.is_3prime_side_of_connection is not False)
             def is_start_3prime(l):
-                return l.get_nt_pos() < 1 and not l.on_fwd_strand
+                is_start = l.get_nt_pos() < 1 and not l.on_fwd_strand
+                return is_start and (l.connection is None or l.is_3prime_side_of_connection is not True)
             def is_end_3prime(l):
-                return l.get_nt_pos() > seg.num_nt-2 and l.on_fwd_strand
+                is_start = l.get_nt_pos() > seg.num_nt-2 and l.on_fwd_strand
+                return is_start and (l.connection is None or l.is_3prime_side_of_connection is not True)
 
             if seg.start5.connection is None:
                 if len(list(filter( is_start_5prime, five_prime_locs ))) == 0:
@@ -3374,6 +3451,9 @@ class SegmentModel(ArbdModel):
                 # if seg.name == "22-2":
                 #     import pdb
                 #     pdb.set_trace()
+                # if _DEBUG_TRACE and seg.name == "69-0" and is_fwd == False:
+                #     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)
                 history.append((seg,pos,end_pos,is_fwd))
@@ -3406,7 +3486,10 @@ class SegmentModel(ArbdModel):
             for l in locs:
                 if _DEBUG_TRACE:
                     print("Tracing 5prime",l)
-                # TODOTODO
+                if l.connection is not None:
+                    print(" Skipping connected 5prime",l)
+                    continue
+                    # TODOTODO
                 pos = seg.contour_to_nt_pos(l.address, round_nt=True)
                 is_fwd = l.on_fwd_strand
                 s = Strand(segname="S{:03d}".format(len(strands)))
@@ -3441,15 +3524,15 @@ class SegmentModel(ArbdModel):
                 return check(last+1)
             else:
                 last = segment.num_nt
-                for sp in segment.strand_pieces[fwd_str]:
-                    if last-sp.end > 1:
+                for sp in reversed(segment.strand_pieces[fwd_str]):
+                    if last-sp.start > 1:
                         return check(last-1)
-                    last = sp.start
+                    last = sp.end
                 return check(last-1)
 
         def add_strand_if_needed(seg,is_fwd):
             history = []
-            if not strands_cover_segment(seg, is_fwd):
+            while not strands_cover_segment(seg, is_fwd):
                 pos = nt = find_nt_not_in_strand(seg, is_fwd)
                 if _DEBUG_TRACE:
                     print("Tracing circular strand", pos, is_fwd)

mrdna/simulate.py

@@ -160,7 +160,7 @@ def multiresolution_simulation( model, output_name,
             model.generate_bead_model( **coarse_model_args )
 
         """ Coarse simulation """
-        simargs = dict(timestep=150e-6 if coarse_local_twist else 200e-6, output_period=coarse_output_period, gpu=gpu)
+        simargs = dict(timestep=100e-6 if coarse_local_twist else 200e-6, output_period=coarse_output_period, gpu=gpu)
         if bond_cutoff > 0 and minimization_steps <= 0:
             run_step(num_steps=0.25*coarse_steps, simargs=simargs)
             model.clear_beads()
@@ -200,7 +200,7 @@ def multiresolution_simulation( model, output_name,
             try:
                 shutil.copytree( get_resource_path("charmm36.nbfix"), "charmm36.nbfix" )
             except FileExistsError:
-                ...
+                pass
 
             model.atomic_simulate( output_name = output_prefix, dry_run = not run_enrg_md )