Commit 6dcc8f64 authored by cmaffeo2's avatar cmaffeo2
Browse files

Made neighborhood for bead model backmapping dynamically increase if needed

parent 72bfa8fb
......@@ -690,20 +690,27 @@ class beadModel():
helixCutoff = dsDnaHelixNeighborDist if b.type[0] == 'd' else ssDnaHelixNeighborDist
allCutoff = dsDnaAllNeighborDist if b.type[0] == 'd' else ssDnaAllNeighborDist
def positionsAreLinear(pos):
if len(pos) == 0: return False
center = np.mean(pos,axis=0)
cPos = pos-center
return np.sum( np.abs( np.linalg.eig( cPos.T.dot(cPos) )[0] ) > 1e-3 ) == 3
ids = []
attempts = 0
while len(ids) <= 3:
posOld = []
while len(ids) <= 3 or not positionsAreLinear(posOld):
if attempts > 15: raise Exception("Too many attempts to find a neighborhood for backmaping bead %d" % b.idx)
ids = simplerModel._getNeighborhoodIds(b, simplerModelCoords, helixCutoff, allCutoff)
posOld = np.array( [simplerModel.particles[i][0].initialPosition for i in ids] )
allCutoff *= 1.2
attempts+=1
posOld = np.array( [simplerModel.particles[i][0].initialPosition for i in ids] )
if attempts > 1:
print("Warning, increased allCutoff to",allCutoff)
posNew = np.array( [simplerModelCoords[i] for i in ids] )
try:
trans[b.idx] = minimizeRmsd( posOld, posNew )
except:
raise Error("Crapola")
# print("ugly")
trans[b.idx] = minimizeRmsd( posOld, posNew )
## Optionally smooth orientations
......
......@@ -762,20 +762,27 @@ class beadModelTwist():
helixCutoff = dsDnaHelixNeighborDist if b.type[0] in ('d','O') else ssDnaHelixNeighborDist
allCutoff = dsDnaAllNeighborDist if b.type[0] in ('d','O') else ssDnaAllNeighborDist
def positionsAreLinear(pos):
if len(pos) == 0: return False
center = np.mean(pos,axis=0)
cPos = pos-center
return np.sum( np.abs( np.linalg.eig( cPos.T.dot(cPos) )[0] ) > 1e-3 ) == 3
ids = []
attempts = 0
while len(ids) <= 3:
posOld = []
while len(ids) <= 3 or not positionsAreLinear(posOld):
if attempts > 15: raise Exception("Too many attempts to find a neighborhood for backmaping bead %d" % b.idx)
ids = simplerModel._getNeighborhoodIds(b, simplerModelCoords, helixCutoff, allCutoff)
posOld = np.array( [simplerModel.particles[i][0].initialPosition for i in ids] )
allCutoff *= 1.2
attempts+=1
posOld = np.array( [simplerModel.particles[i][0].initialPosition for i in ids] )
if attempts > 1:
print("Warning, increased allCutoff to",allCutoff)
posNew = np.array( [simplerModelCoords[i] for i in ids] )
try:
trans[b.idx] = minimizeRmsd( posOld, posNew )
except:
raise Exception("Failed to find orientation of atom %d in the coarser model" % b.idx)
# print("ugly")
trans[b.idx] = minimizeRmsd( posOld, posNew )
## Optionally smooth orientations
......
Supports Markdown
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