Source code for arbdmodel.coords

from .logger import logger
import numpy as np

def minimizeRmsd(coordsB, coordsA, weights=None, maxIter=100):
    ## Going through many iterations wasn't really needed
    tol = 1
    count = 0

    R = np.eye(3)
    comB = np.zeros([3,])
    cNext = coordsB

    while tol > 1e-6:
        q,cB,comA = _minimizeRmsd(cNext,coordsA, weights)
        R = R.dot(quaternion_to_matrix(q))
        assert( np.all(np.isreal( R )) )

        comB += cB
        cLast = cNext
        cNext = (coordsB-comB).dot(R)

        tol = np.sum(((cNext-cLast)**2)[:]) / np.max(np.shape(coordsB))
        if count > maxIter:
            Exception("Exceeded maxIter (%d)" % maxIter)
        count += 1

    print("%d iterations",count)
    return R, comB, comA


[docs] def minimizeRmsd(coordsB, coordsA, weights=None): q,comB,comA = _minimizeRmsd(coordsB, coordsA, weights) assert( np.all(np.isreal( q )) ) return quaternion_to_matrix(q),comB,comA
## http://onlinelibrary.wiley.com/doi/10.1002/jcc.21439/full def _minimizeRmsd(coordsB, coordsA, weights=None): A = coordsA B = coordsB shapeA,shapeB = [np.shape(X) for X in (A,B)] for s in (shapeA,shapeB): assert( len(s) == 2 ) A,B = [X.T if s[1] > s[0] else X for X,s in zip([A,B],(shapeA,shapeB))] # TODO: print warning shapeA,shapeB = [np.shape(X) for X in (A,B)] assert( shapeA == shapeB ) for X,s in zip((A,B),(shapeA,shapeB)): assert( s[1] == 3 and s[0] >= s[1] ) # if weights is None: weights = np.ones(len(A)) if weights is None: comA,comB = [np.mean( X, axis=0 ) for X in (A,B)] else: assert( len(weights[:]) == len(B) ) W = np.diag(weights) comA,comB = [np.sum( W.dot(X), axis=0 ) / np.sum(W) for X in (A,B)] A = np.array( A-comA ) B = np.array( B-comB ) if weights is None: s = A.T.dot(B) else: s = A.T.dot(W.dot(B)) sxx,sxy,sxz = s[0,:] syx,syy,syz = s[1,:] szx,szy,szz = s[2,:] K = [[ sxx+syy+szz, syz-szy, szx-sxz, sxy-syx], [syz-szy, sxx-syy-szz, sxy+syx, sxz+szx], [szx-sxz, sxy+syx, -sxx+syy-szz, syz+szy], [sxy-syx, sxz+szx, syz+szy, -sxx-syy+szz]] K = np.array(K) # GA = np.trace( A.T.dot(W.dot(A)) ) # GB = np.trace( B.T.dot(W.dot(B)) ) ## Finding GA/GB can be done more quickly # I = np.eye(4) # x0 = (GA+GB)*0.5 # vals = newtoon(lambda x: np.det(K-x*I), x0 = x0) vals, vecs = np.linalg.eig(K) i = np.argmax(vals) q = vecs[:,i] # RMSD = np.sqrt( (GA+GB-2*vals[i]) / len(A) ) # print("CHECK:", K.dot(q)-vals[i]*q ) return q, comB, comA
[docs] def quaternion_to_matrix(q): assert(len(q) == 4) ## It looks like the wikipedia article I used employed a less common convention for q (see below ## http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions#Rotation_matrix_.E2.86.94_quaternion # q1,q2,q3,q4 = q # R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q4), 2*(q1*q3 + q2*q4)], # [ 2*(q1*q2 + q3*q4), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q4)], # [ 2*(q1*q3 - q2*q4), 2*(q1*q4 + q2*q3), 1-2*(q2*q2 + q1*q1)]] q = q / np.linalg.norm(q) q0,q1,q2,q3 = q R = [[1-2*(q2*q2 + q3*q3), 2*(q1*q2 - q3*q0), 2*(q1*q3 + q2*q0)], [ 2*(q1*q2 + q3*q0), 1-2*(q1*q1 + q3*q3), 2*(q2*q3 - q1*q0)], [ 2*(q1*q3 - q2*q0), 2*(q1*q0 + q2*q3), 1-2*(q2*q2 + q1*q1)]] return np.array(R)
[docs] def quaternion_from_matrix( R ): R=R.T q = np.empty(4) if R[2,2] < 0: if R[0,0] > R[1,1]: trace = 1.0 + R[0,0] - R[1,1] - R[2,2] s = 2.0 * np.sqrt(trace) if R[1,2] < R[2,1]: s = -s q[0] = (R[1,2] - R[2,1]) / s q[1] = 0.25 * s q[2] = (R[0,1] + R[1,0]) / s q[3] = (R[2,0] + R[0,2]) / s if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 1],0)): q[1] = 1 else: trace = 1.0 - R[0,0] + R[1,1] - R[2,2] s = 2.0 * np.sqrt(trace) if R[2,0] < R[0,2]: s = -s q[0] = (R[2,0] - R[0,2]) / s q[1] = (R[0,1] + R[1,0]) / s q[2] = 0.25 * s q[3] = (R[1,2] + R[2,1]) / s if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 2],0)): q[2] = 1 else: if R[0,0] < -R[1,1]: trace = 1.0 - R[0,0] - R[1,1] + R[2,2] s = 2.0 * np.sqrt(trace) if R[0,1] < R[1,0]: s = -s q[0] = (R[0,1] - R[1,0]) / s q[1] = (R[2,0] + R[0,2]) / s q[2] = (R[1,2] + R[2,1]) / s q[3] = 0.25 * s if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 3],0)): q[3] = 1 else: trace = 1.0 + R[0,0] + R[1,1] + R[2,2] s = 2.0 * np.sqrt(trace) q[0] = 0.25 * s q[1] = (R[1,2] - R[2,1]) / s q[2] = (R[2,0] - R[0,2]) / s q[3] = (R[0,1] - R[1,0]) / s if np.isclose(trace,1) and np.all(np.isclose([x for i,x in enumerate(q) if i != 0],0)): q[0] = 1 assert( q[0] >= 0 ) return q
def __quaternion_from_matrix__deprecated( R ): e1 = R[0] e2 = R[1] e3 = R[2] # d1 = 0.5 * np.sqrt( 1+R[0,0]+R[1,1]+R[2,2] ) # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] ) # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] ) # d2 = 0.5 * np.sqrt( 1+R[0,0]-R[1,1]-R[2,2] ) d1 = 1+R[0,0]+R[1,1]+R[2,2] d2 = 1+R[0,0]-R[1,1]-R[2,2] d3 = 1-R[0,0]+R[1,1]-R[2,2] d4 = 1-R[0,0]-R[1,1]+R[2,2] maxD = max((d1,d2,d3,d4)) d = 0.5 / np.sqrt(maxD) if d1 == maxD: return np.array(( 1.0/(4*d), d * (R[2,1]-R[1,2]), d * (R[0,2]-R[2,0]), d * (R[1,0]-R[0,1]) )) elif d2 == maxD: return np.array(( d * (R[2,1]-R[1,2]), 1.0/(4*d), d * (R[0,1]+R[1,0]), d * (R[0,2]+R[2,0]) )) elif d3 == maxD: return np.array(( d * (R[0,2]-R[2,0]), d * (R[0,1]+R[1,0]), 1.0/(4*d), d * (R[1,2]+R[2,1]) )) elif d4 == maxD: return np.array(( d * (R[1,0]-R[0,1]), d * (R[0,2]+R[2,0]), d * (R[1,2]+R[2,1]), 1.0/(4*d) ))
[docs] def quaternion_product(a, b): assert(len(a) == 4) assert(len(b) == 4) ab = np.empty(4) ab[0] = a[0]*b[0]-a[1:].dot(b[1:]) ab[1:] = a[0]*b[1:] + b[0]*a[1:] + np.cross(a[1:],b[1:]) return ab
[docs] def quaternion_inverse(q): assert(len(q) == 4) qinv = np.array(q) qinv[1:] = -qinv[1:] qinv = qinv / (q[0]**2 + q[1:].dot(q[1:])) return qinv
[docs] def quaternion_exp(q,t): assert(len(q) == 4) omega = np.arccos(q[0]) v = q[1:] v = v/np.linalg.norm(v) qexp = np.empty(4) qexp[0] = np.cos(omega*t) qexp[1:] = np.sin(omega*t)*v return qexp
[docs] def quaternion_slerp(q1,q2,t): assert(len(q1) == 4) assert(len(q2) == 4) assert(t >= 0 and t <= 1) q1_inv = quaternion_inverse(q1) return quaternion_product( q1, quaternion_exp( quaternion_product( q1_inv, q2 ), t ) )
[docs] def rotationAboutAxis(axis,angle, normalizeAxis=True): if normalizeAxis: axis = axis / np.linalg.norm(axis) angle = angle * 0.5 * np.pi/180 cos = np.cos( angle ) sin = np.sin( angle ) q = [cos] + [sin*x for x in axis] return quaternion_to_matrix(q)
# By Chun
[docs] def Generate_spanning_vectors(bv1, bv2, bv3, dimensions, buff=5): dd = max(dimensions) + 2 * buff n1 = round(np.linalg.norm(bv1)/dd) + 1 n2 = round(np.linalg.norm(bv2)/dd) + 1 n3 = round(np.linalg.norm(bv3)/dd) + 1 v1 = np.array(bv1) /n1 v2 = np.array(bv2) /n2 v3 = np.array(bv3) /n3 return v1, v2, v3, n1, n2, n3
[docs] def Generate_coordinates(bv1, bv2, bv3, n1, n2, n3, num_copy, origin, replica_index): ori_vec = np.array(origin) if n1 * n2 * n3 > num_copy: check = True else: check = False np.random.seed(42 + replica_index) count = 0 inds = [] while count < num_copy: ind = (np.random.randint(0, n1), np.random.randint(0, n2), np.random.randint(0, n3)) if not check: inds.append(ind) count += 1 elif check: if not ind in inds: inds.append(ind) count += 1 else: while ind in inds: ind = (np.random.randint(0, n1), np.random.randint(0, n2), np.random.randint(0, n3)) inds.append(ind) count += 1 coors = [] for ind in inds: r1 = (ind[0] - (n1 - 1) * 0.5) * bv1 r2 = (ind[1] - (n2 - 1) * 0.5) * bv2 r3 = (ind[2] - (n3 - 1) * 0.5) * bv3 coors.append((r1 + r2 + r3 + ori_vec).tolist()) return coors
#Code to here from simplearbd
[docs] def readArbdCoords(fname): logger.warning('readArbdCoords is deprecated. Please update your code to use read_arbd_coordinates') return read_arbd_coordinates(fname)
[docs] def read_arbd_coordinates(fname): coords = [] with open(fname) as fh: for line in fh: coords.append([float(x) for x in line.split()[1:]]) return np.array(coords)
[docs] def readAvgArbdCoords(*args,**kwargs): logger.warning('readAvgArbdCoords is deprecated. Please update your code to use read_average_arbd_coordinates') return read_average_arbd_coordinates(*args,**kwargs)
[docs] def read_average_arbd_coordinates(psf,pdb,dcd,rmsd_threshold=3.5, first_frame=0, last_frame=-1, stride=1): import MDAnalysis as mda usel = mda.Universe(psf, dcd) sel = usel.select_atoms("name D*") # r0 = ref.xyz[0,ids,:] ts = usel.trajectory[last_frame] r0 = sel.positions pos = [] for t in range(ts.frame-1,first_frame-1,-stride): usel.trajectory[t] R,comA,comB = minimizeRmsd(sel.positions,r0) r = np.array( [(r-comA).dot(R)+comB for r in sel.positions] ) rmsd = np.mean( (r-r0)**2 ) r = np.array( [(r-comA).dot(R)+comB for r in usel.atoms.positions] ) pos.append( r ) if rmsd > rmsd_threshold**2: break t0=t+1 logger.info(f"Averaging coordinates in {dcd} after frame {t0}") pos = np.mean(pos, axis=0) return pos
[docs] def calculate_dimensions_from_cell_vectors(cell_vectors, cell_origin=None, buffer_factor=1.2): """ Calculate simulation box dimensions from cell basis vectors with buffer. Args: cell_vectors: List of 3 cell basis vectors [[x1,y1,z1], [x2,y2,z2], [x3,y3,z3]] cell_origin: Cell origin coordinates [x,y,z], defaults to [0,0,0] buffer_factor: Factor to scale the dimensions for additional buffer space Returns: dimensions: 3D box dimensions as a tuple (x_dim, y_dim, z_dim) """ if cell_vectors is None: return None # Make sure we have proper cell vectors if len(cell_vectors) != 3: raise ValueError("Expected 3 cell basis vectors, got {}".format(len(cell_vectors))) # Convert to numpy arrays for easier manipulation vectors = [np.array(v) for v in cell_vectors] # Calculate maximum extent along each dimension dimensions = [0, 0, 0] for i in range(3): for v in vectors: dimensions[i] += abs(v[i]) # Apply buffer factor dimensions = [dim * buffer_factor for dim in dimensions] return tuple(dimensions)
[docs] def unit_quat_conversions(): for axis in [[0,0,1],[1,1,1],[1,0,0],[-1,-2,0]]: for angle in np.linspace(-180,180,10): R = rotationAboutAxis(axis, angle) R2 = quaternion_to_matrix( quaternion_from_matrix( R ) ) if not np.all( np.abs(R-R2) < 0.01 ): import pdb pdb.set_trace() quaternion_to_matrix( quaternion_from_matrix( R ) )
if __name__ == "__main__": unit_quat_conversions()