Commit 671f1605 authored by cmaffeo2's avatar cmaffeo2
Browse files

Fixed multiple issues; 3-prime ends of nucleotides now terminate correctly;...

Fixed multiple issues; 3-prime ends of nucleotides now terminate correctly; strands not cleared unless atoms exist; strands regenerated with atomic model; cadnano models initialized with more correct azimuthal angle
parent d9626894
......@@ -1105,13 +1105,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.pdb
# fixedAtomsCol B
# minimize 2400
# fixedAtoms off
minimize 2400
}} else {{
bincoordinates {output_directory}/$prefix-$nLast.restart.coor
......
......@@ -3,8 +3,8 @@ import scipy.optimize as opt
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter as savgol
from .nonbonded import NonbondedScheme
from .. import get_resource_path
from mrdna.model.nonbonded import NonbondedScheme
from mrdna import get_resource_path
dxParam = -1
......@@ -33,7 +33,7 @@ y0 = d[:,0] # kcal/mol for 10bp with 10bp, according to JY
y0 = y0 - np.mean(y0[x0>(38.5 + dxParam)]) # set 0 to tail of potential
y0 = savgol(y0, int(5.0/(x0[1]-x0[0])), 3) # smooth potential
def pot(x,*parms):
def parametric_potential(x,*parms):
x = np.array(x)
u = interp1d(x0, np.array(parms),
bounds_error = False, fill_value="extrapolate",
......@@ -41,7 +41,7 @@ def pot(x,*parms):
u[x>np.max(x0)] = 0
return u
def potForBp(potfn,x,numBp,*parms):
def integrated_potential(potfn,x,numBp,*parms):
u = np.zeros(np.shape(x))
for h in range(-numBp,numBp+1):
r = np.sqrt(x**2 + (3.4*h)**2)
......@@ -50,7 +50,8 @@ def potForBp(potfn,x,numBp,*parms):
return u
def fitFun(x,*parms):
return 10 * potForBp(pot,x,20,*parms)
## 1-turn iteracting with 4 turns (more is not needed)
return 10 * integrated_potential(parametric_potential,x,20,*parms)
xinterp = np.linspace(np.min(x0),np.max(x0),1000) # interpolate potentialb
yinterp = interp1d(x0,y0)(xinterp)
......@@ -62,9 +63,9 @@ def nbPot(x,bps1,bps2):
if larger >= 2:
largerInt = int(0.5*larger)
smaller *= larger/(2*largerInt)
y = smaller * potForBp(pot,x,largerInt,*popt)
y = smaller * integrated_potential(parametric_potential,x,largerInt,*popt)
else:
y = smaller * larger * pot(x,*popt)
y = smaller * larger * parametric_potential(x,*popt)
y = maxForce(x,y, 100)
y = savgol(y, int(5.0/(x0[1]-x0[0])), 3)
return y
......@@ -99,3 +100,10 @@ class nbDnaScheme(NonbondedScheme):
u = nbPot(r, 0.5*typeA.nts, 0.5*typeB.nts)
return u
nbDnaScheme = nbDnaScheme()
if __name__ == "__main__":
## run some tests
np.savetxt("jj-filtered.dat", np.array((x0,y0)).T)
y = fitFun(x,*popt)
np.savetxt("1-bp-fit.dat", np.array((x,y)).T),
......@@ -5,7 +5,7 @@ import os,sys
from glob import glob
import re
from ..coords import readArbdCoords, readAvgArbdCoords
from ..coords import readArbdCoords, readAvgArbdCoords, rotationAboutAxis
from ..segmentmodel import SegmentModel, SingleStrandedSegment, DoubleStrandedSegment
from ..model.dna_sequence import m13 as m13seq
......@@ -311,7 +311,8 @@ class cadnano_part(SegmentModel):
## TODO get sequence from cadnano api
if zMid in strandOccupancies[0] and zMid in strandOccupancies[1]:
kwargs['num_bp'] = numBps
seg = DoubleStrandedSegment(**kwargs,**posargs1)
start_orientation = rotationAboutAxis(np.array((0,0,1)), (zid1-3.25)*34)
seg = DoubleStrandedSegment(**kwargs,**posargs1, start_orientation = start_orientation)
elif zMid in strandOccupancies[0]:
kwargs['num_nt'] = numBps
seg = SingleStrandedSegment(**kwargs,**posargs1)
......
......@@ -506,6 +506,9 @@ class Segment(ConnectableElement, Group):
else:
orientation0 = np.eye(3) if axis.dot(zAxis) > 0 else \
rotationAboutAxis( np.array((1,0,0)), 180, normalizeAxis=False )
if self.start_orientation is not None:
orientation0 = orientation0.dot(self.start_orientation)
orientation = rotationAboutAxis( axis, self.twist_per_nt*self.contour_to_nt_pos(s), normalizeAxis=False )
orientation = orientation.dot(orientation0)
else:
......@@ -572,7 +575,7 @@ class Segment(ConnectableElement, Group):
if scale is not None and scale != 1:
for a in atoms:
a.position = scale*a.position
a.beta = 0
a.fixed = 1
atoms.position = pos - atoms.atoms_by_name["C1'"].collapsedPosition()
else:
if scale is not None and scale != 1:
......@@ -583,7 +586,7 @@ class Segment(ConnectableElement, Group):
for a in atoms:
if a.name[-1] in ("'","P","T"):
a.position = scale*(a.position-r0) + r0
a.beta = 0
a.fixed = 1
atoms.position = pos
atoms.contour_position = contour_position
......@@ -1321,14 +1324,15 @@ class Strand(Group):
# pdb.set_trace()
# assert(s.end != s.start)
assert( s.num_nt == 1 or (np.linalg.norm( seg.contour_to_position(contour[-1]) - seg.contour_to_position(contour[0]) ) > 0.1) )
nucleotide_count = 0
for c,seq in zip(contour,s.get_sequence()):
nucleotide_count += 1
if last is None and not self.is_circular:
seq = "5"+seq
if strand_segment_count == len(s.strand_segments) and c == 1 and not self.is_circular:
if strand_segment_count == len(s.strand_segments) and nucleotide_count == s.num_nt and not self.is_circular:
seq = seq+"3"
nt = seg._generate_atomic_nucleotide( c, s.is_fwd, seq, scale, s )
s.add(nt)
## Join last basepairs
if last is not None:
......@@ -1515,6 +1519,7 @@ class SegmentModel(ArbdModel):
for seg in self.segments:
for d in ('fwd','rev'):
seg.strand_pieces[d] = []
self._generate_strands()
def clear_beads(self):
return self._clear_beads()
......@@ -1525,7 +1530,8 @@ class SegmentModel(ArbdModel):
s.clear_all()
self.clear_all(keep_children=True)
self.clear_atomic()
if len(self.strands[0].children[0].children) > 0:
self.clear_atomic()
## Check that it worked
assert( len([b for b in self]) == 0 )
......@@ -2355,6 +2361,9 @@ class SegmentModel(ArbdModel):
try:
for s in self.strands:
self.children.remove(s)
for seg in self.segments:
for d in ('fwd','rev'):
seg.strand_pieces[d] = []
except:
...
self.strands = strands = []
......@@ -2674,7 +2683,7 @@ class SegmentModel(ArbdModel):
# print("WARNING: could not find 'P' atom in {}:{} or {}:{}".format( segI, ntI_idx, segJ, ntJ_idx ))
...
print("PUSH BONDS:", len(push_bonds))
# print("PUSH BONDS:", len(push_bonds))
if not self.useTclForces:
with open("%s.exb" % output_name, 'a') as fh:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment