|
| 1 | +import numpy as np |
| 2 | + |
| 3 | +from scipy.interpolate import InterpolatedUnivariateSpline |
| 4 | + |
| 5 | + |
| 6 | +def interpolate_chain(positions, extra_bonds=[], target_N=90000, fine_grain=10): |
| 7 | + """ |
| 8 | + Converts an individual polymer of any length to a smoothed chain with fixed distance between neighboring monomers. |
| 9 | + Proceeds by cubic spline interpolation as follows. |
| 10 | +
|
| 11 | + 1. Interpolate the data using cubic spline |
| 12 | + 2. Evaluate cubic spline at targetN*10 values |
| 13 | + 3. Rescale the evaluated spline such that total distance is targetN |
| 14 | + 4. Select targetN points along the path with distance between neighboring points _along the chain_ equal to 1. |
| 15 | + 5. Select appropriate bonds along the interpolated path to mimic the original chain topology |
| 16 | +
|
| 17 | + Parameters |
| 18 | + ---------- |
| 19 | + positions : Nx3 float array |
| 20 | + Input xyz coordinates |
| 21 | + extra_bonds : Mx2 int array |
| 22 | + Non-contiguous extra bonds to be interpolated (optional) |
| 23 | + target_N : int |
| 24 | + Length of output polymer. |
| 25 | + It is not advised to make it many times less than N |
| 26 | +
|
| 27 | + Returns |
| 28 | + ------- |
| 29 | + (about target_N) x 3 float array of interpolated positions |
| 30 | + (about target_N-1 + M) x 2 int array of interpolated bonds |
| 31 | + """ |
| 32 | + |
| 33 | + N = len(positions) |
| 34 | + target_data_size = target_N*fine_grain |
| 35 | + |
| 36 | + eval_range = np.arange(N) |
| 37 | + target_range = np.arange(0, N-1, N/float(target_data_size)) |
| 38 | + |
| 39 | + splined = np.zeros((len(target_range), 3), float) |
| 40 | + |
| 41 | + for i in range(3): |
| 42 | + spline = InterpolatedUnivariateSpline(eval_range, positions[:, i], k=3) |
| 43 | + splined[:, i] = spline(target_range) |
| 44 | + |
| 45 | + pos = np.arange(1, target_N) |
| 46 | + dists = np.linalg.norm(np.diff(splined, axis=0), axis=1) |
| 47 | + |
| 48 | + cum_dists = np.cumsum(dists) |
| 49 | + |
| 50 | + splined /= (cum_dists[-1] / N) |
| 51 | + cum_dists /= (cum_dists[-1] / target_N) |
| 52 | + |
| 53 | + searched_pos = np.searchsorted(cum_dists, pos) |
| 54 | + searched_bonds = np.searchsorted(target_range[searched_pos], extra_bonds) |
| 55 | + |
| 56 | + v1 = cum_dists[searched_pos] |
| 57 | + v2 = cum_dists[searched_pos-1] |
| 58 | + |
| 59 | + p1 = (v1-np.floor(v1)) / (v1-v2) |
| 60 | + p2 = 1-p1 |
| 61 | + |
| 62 | + interp_pos = p2[:, None]*splined[searched_pos] + p1[:, None]*splined[searched_pos-1] |
| 63 | + interp_bonds = np.asarray(list(zip(pos[:-1]-1, pos[:-1])) + searched_bonds.tolist()) |
| 64 | + |
| 65 | + return interp_pos, interp_bonds |
0 commit comments