Reinstated fixed atoms for atomic simulations

mrdna/model/arbdmodel.py

@@ -481,7 +481,7 @@ class PdbModel(Transformable, Parent):
     def _updateParticleOrder(self):
         pass
 
-    def writePdb(self, filename):
+    def writePdb(self, filename, beta_from_fixed=False):
         if self.cacheInvalid:
             self._updateParticleOrder()
         with open(filename,'w') as fh:
@@ -501,6 +501,9 @@ class PdbModel(Transformable, Parent):
                     idx = "{:>6d}".format(idx)
                 data['idx'] = idx
 
+                if beta_from_fixed:
+                    data['beta'] = 1 if 'fixed' in p.__dict__ else 0
+
                 pos = p.collapsedPosition()
                 dig = [max(int(np.log10(np.abs(x)+1e-6)//1),0)+1 for x in pos]
                 for d in dig: assert( d <= 7 )
@@ -1105,12 +1108,12 @@ cellBasisVector3 0 0 {cellZ}
 
 if {{$nLast == 0}} {{
     temperature 300
-    # fixedAtoms on
-    # fixedAtomsForces on
-    # fixedAtomsFile $prefix.pdb
-    # fixedAtomsCol B
-    # minimize 2400
-    # fixedAtoms off
+    fixedAtoms on
+    fixedAtomsForces on
+    fixedAtomsFile $prefix.fixed.pdb
+    fixedAtomsCol B
+    minimize 2400
+    fixedAtoms off
     minimize 2400
 }} else {{
     bincoordinates  {output_directory}/$prefix-$nLast.restart.coor
@@ -1127,6 +1130,7 @@ run {num_steps:d}
 
         if output_directory == '': output_directory='.'
         self.writePdb( output_name + ".pdb" )
+        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 )
         os.sync()

mrdna/readers/segmentmodel_from_lists.py

@@ -70,18 +70,6 @@ def basepairs_and_stacks_to_helixmap(basepairs,stacks_above):
     return helixmap, helixrank, is_fwd
 
 
-def SegmentModelFromPdb(*args,**kwargs):
-    u = mda.Universe(*args,**kwargs)
-    for a in u.select_atoms("resname DT* and name C7").atoms:
-        a.name = "C5M"
-
-    ## Find basepairs and base stacks
-    centers,transforms = find_base_position_orientation(u)
-    bps = find_basepairs(u, centers, transforms)
-    stacks = find_stacks(u, centers, transforms)
-
-    ## Build map from residue index to helix index
-
 def set_splines(seg, coordinates, hid, hmap, hrank, fwd, orientation=None):
     maxrank = np.max( hrank[hmap==hid] )
     if maxrank == 0:

mrdna/readers/segmentmodel_from_pdb.py

@@ -258,6 +258,8 @@ def SegmentModelFromPdb(*args,**kwargs):
     bps = find_basepairs(u, centers, transforms)
     stacks = find_stacks(u, centers, transforms)
 
+    pdb.set_trace()
+
     ## Find three-prime ends
     three_prime = -np.ones(bps.shape, dtype=np.int)
     for s in u.segments:

mrdna/segmentmodel.py

@@ -638,7 +638,6 @@ class Segment(ConnectableElement, Group):
             if scale is not None and scale != 1:
                 for a in atoms:
                     a.position = scale*a.position
-                    a.fixed = 1
             atoms.position = pos - atoms.atoms_by_name["C1'"].collapsedPosition()
         else:
             if scale is not None and scale != 1:
@@ -649,6 +648,7 @@ class Segment(ConnectableElement, Group):
                 for a in atoms:
                     if a.name[-1] in ("'","P","T"):
                         a.position = scale*(a.position-r0) + r0
+                    else:
                         a.fixed = 1
             atoms.position = pos