Skip to content

Commit f283f31

Browse files
anyuzxanyuzx
authored
Merge pull request #5 from anyuzx/dynamics
Add dynamic prediction functionality
2 parents ea97db8 + 9439bfb commit f283f31

File tree

1 file changed

+238
-3
lines changed

1 file changed

+238
-3
lines changed

‎HippsDimes.py

Lines changed: 238 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,127 @@
3333

3434
#------------------------------------------------------------------#
3535
# Helper functions
36-
36+
def compute_acf_general_theory(i, j, t, a, zeta=1.0):
37+
"""
38+
Numerically compute the autocorrelation function (ACF) for monomers i, j
39+
using the connectivity matrix `a`. Returns both the time-dependent
40+
ACF and the corresponding MSD for each time point.
41+
42+
Parameters
43+
----------
44+
i, j : int
45+
Indices of the monomers for which to calculate the correlation function.
46+
t : array_like
47+
A 1D array of time points (lag times).
48+
a : np.ndarray
49+
The connectivity (or "Laplacian") matrix for the polymer/chain.
50+
zeta : float
51+
Friction coefficient (if part of the model). Default is 1.0.
52+
53+
Returns
54+
-------
55+
output : np.ndarray
56+
2D array. First column is time `t`, second column is the ACF values at each time.
57+
msd : np.ndarray
58+
1D array of the same length as `t`, giving the mean-square displacement
59+
inferred from the ACF. This is returned as a second column alongside `t`.
60+
"""
61+
eigvalue, eigvector = scipy.linalg.eigh(a)
62+
eigvalue_inv = 1.0 / eigvalue
63+
64+
# difference in eigenvector components for monomers i and j
65+
vpi_vpj = eigvector[i, :] - eigvector[j, :]
66+
67+
# normal_modes_square_mean = -(1 / eigenvalue) but filter out any inf
68+
normal_modes_square_mean = - np.nan_to_num(eigvalue_inv, posinf=0.0, neginf=0.0)
69+
70+
# Expand time dimension for broadcast
71+
t_reshaped = np.expand_dims(t, axis=-1)
72+
# Effective relaxation times
73+
tau_p = - zeta / eigvalue
74+
decay_factor = np.exp(-t_reshaped / tau_p)
75+
76+
# ACF(t)
77+
res = 3.0 * np.sum(vpi_vpj**2 * decay_factor * normal_modes_square_mean, axis=-1)
78+
# Equilibrium part
79+
res_eq = 3.0 * np.sum(vpi_vpj**2 * normal_modes_square_mean, axis=-1)
80+
81+
# Combine results: first column is time, second is res
82+
two_point_acf = np.column_stack((t, res))
83+
84+
# The mean-square displacement from the correlation function
85+
two_point_msd = 2.0 * (res_eq - two_point_acf[:, 1])
86+
two_point_msd = np.column_stack((t, two_point_msd))
87+
88+
return two_point_acf, two_point_msd
89+
90+
def compute_m1_general_theory(i, t, a, zeta=1.0):
91+
"""
92+
Compute the single-monomer mean-square displacement (MSD) for monomer i,
93+
given the connectivity matrix `a`.
94+
95+
Parameters
96+
----------
97+
i : int
98+
Index of the monomer.
99+
t : array_like
100+
A 1D array of time points (lag times).
101+
a : np.ndarray
102+
The connectivity (or "Laplacian") matrix for the polymer/chain.
103+
zeta : float
104+
Friction coefficient. Default is 1.0.
105+
106+
Returns
107+
-------
108+
msd : np.ndarray
109+
2D array. First column is time `t`, second column is the MSD for
110+
monomer i at those times.
111+
"""
112+
eigvalue, eigvector = scipy.linalg.eigh(a)
113+
eigvalue_inv = 1.0 / eigvalue
114+
vpi = eigvector[i, :]
115+
116+
# Filter out infinities
117+
normal_modes_square_mean = - np.nan_to_num(eigvalue_inv, posinf=0.0, neginf=0.0)
118+
119+
# Expand time dimension for broadcast
120+
t_reshaped = np.expand_dims(t, axis=-1)
121+
tau_p = - zeta / eigvalue
122+
decay_factor = np.exp(-t_reshaped / tau_p)
123+
124+
# The time-dependent part
125+
res = 3.0 * np.sum(vpi**2 * decay_factor * normal_modes_square_mean, axis=-1)
126+
# Equilibrium radius
127+
r2_eq = 3.0 * np.sum(vpi**2 * normal_modes_square_mean, axis=-1)
128+
# MSD
129+
msd_data = 2.0 * (r2_eq - res)
130+
131+
# Combine time with MSD
132+
msd = np.column_stack((t, msd_data))
133+
return msd
134+
135+
def Ornstein_Uhlenbeck_update(x, dt, k, zeta, beta, b = 0.0, method='euler-maruyama'):
136+
"""
137+
Update variable x for a Ornstein Uhlenbeck process
138+
x: Array for value of x of each degree of freedom
139+
k: Array for spring constant for each degree of freedom
140+
zeta: one value
141+
beta: one value
142+
"""
143+
if isinstance(x, np.ndarray):
144+
rand_noise = np.random.randn(*x.shape)
145+
else:
146+
rand_noise = np.random.randn()
147+
148+
if method == 'euler-maruyama':
149+
dx = - k[:, np.newaxis] * x * dt / zeta + b * dt / zeta + np.sqrt(2.0 * dt / (zeta * beta)) * rand_noise
150+
x_new = x + dx
151+
elif method == 'exact':
152+
theta = k[:, np.newaxis] / zeta
153+
sigma = (2. / (zeta * beta)) ** .5
154+
mu = np.exp(- theta * dt)
155+
x_new = x * mu + np.nan_to_num(np.sqrt((sigma ** 2. / (2. * theta)) * (1. - mu ** 2.))) * rand_noise
156+
return x_new
37157

38158
def construct_connectivity_matrix_rouse(n, k):
39159
"""
@@ -244,8 +364,6 @@ def objective_func(rc, A_mtx, cmap_exp):
244364
return res
245365

246366
# FUNCTION TO CONVERT CMAP TO DMAP
247-
248-
249367
def cmap2dmap_core(cmap_exp, rc, alpha, not_normalize, norm_max=1.0, mode='log'):
250368
# rc is the prefactor
251369
# norm_max is the maximum contact probability
@@ -445,6 +563,123 @@ def run(self, epoch, general_method='optimization', **kwargs):
445563

446564
return loss_array, dmap_maxent, self.A
447565

566+
class Dynamics:
567+
def __init__(self, input, M=None, k=None, model=None):
568+
if isinstance(input, int) and M is None and k is not None:
569+
if not isinstance(k, int) and not isinstance(k, float):
570+
sys.stdout.write('Spring constant should be a number')
571+
sys.exit(0)
572+
if model != 'rouse' and model is not None:
573+
sys.stdout.write("Please specify model to be 'rouse'")
574+
sys.exit(0)
575+
576+
self.A = construct_connectivity_matrix_rouse(input, k)
577+
self.eigvalue, self.eigvector = scipy.linalg.eigh(self.A)
578+
self.N = input
579+
elif isinstance(input, int) and M is not None and k is not None:
580+
if not isinstance(k, int) and not isinstance(k, float):
581+
sys.stdout.write('Spring constant should be a number')
582+
sys.exit(0)
583+
if not isinstance(M, int):
584+
sys.stdout.write('Number of random cross links need to be an integer')
585+
sys.exit(0)
586+
if model != 'random':
587+
sys.stdout.write("Please specify model to be 'random'")
588+
sys.exit(0)
589+
self.A = construct_connectivity_matrix_random(input, M, k)
590+
self.eigvalue, self.eigvector = scipy.linalg.eigh(self.A)
591+
self.N = input
592+
elif isinstance(input, np.ndarray) and M is None and k is None:
593+
if len(input.shape) !=2 or input.shape[0] != input.shape[1]:
594+
sys.stdout.write('The connectivity matrix should be a square matrix')
595+
sys.exit(0)
596+
if not np.allclose(input, input.T):
597+
sys.stdout.write('The connectivity matrix should be a symmetrix real matrix')
598+
sys.exit(0)
599+
600+
self.A = input
601+
self.eigvalue, self.eigvector = scipy.linalg.eigh(self.A)
602+
self.N = input.shape[0]
603+
604+
def generateXYZ(self, force_positive_definite = False):
605+
self.xyz = a2xyz(self.A, force_positive_definite = force_positive_definite)
606+
self.modes = self.eigvector.T @ self.xyz
607+
608+
def initialize(self, dt, zeta, beta):
609+
if not isinstance(dt, int) and not isinstance(dt, float):
610+
sys.stdout.write('Time step should be a number')
611+
sys.exit(0)
612+
if not isinstance(zeta, int) and not isinstance(zeta, float):
613+
sys.stdout.write('Friction coefficient step should be a number')
614+
sys.exit(0)
615+
if not isinstance(beta, int) and not isinstance(beta, float):
616+
sys.stdout.write('Temperature step should be a number')
617+
sys.exit(0)
618+
elif beta <= 0.0:
619+
sys.stdout.write('Temperature should be positive')
620+
sys.exit(0)
621+
622+
self.zeta = zeta
623+
self.beta = beta
624+
self.dt = dt
625+
626+
def updateModes(self, method='euler-maruyama'):
627+
try:
628+
self.zeta
629+
self.beta
630+
self.dt
631+
except AttributeError:
632+
sys.stdout.write('Please run initialize() first')
633+
sys.exit(0)
634+
635+
self.modes = Ornstein_Uhlenbeck_update(self.modes, self.dt, - self.eigvalue, self.zeta, self.beta, method=method)
636+
# self.modes = OU.OU(self.modes, self.dt, - self.eigvalue, self.zeta, self.beta)
637+
638+
def updateXYZ(self):
639+
self.xyz = self.eigvector @ self.modes
640+
641+
def run(self, T, update=1, every=1, initial_conformation=None, method='euler-maruyama'):
642+
"""
643+
T: number of timesteps
644+
update: update x,y,z positions every this many timesteps
645+
every: save x,y,z positions to the trajectory every this many timesteps
646+
initial_conformation: initial conformation of the simulation
647+
"""
648+
if not isinstance(T, int):
649+
sys.stdout.write('Number of steps should be an integer')
650+
sys.exit(0)
651+
652+
if initial_conformation is None:
653+
try:
654+
self.xyz
655+
self.modes
656+
except AttributeError:
657+
self.generateXYZ()
658+
else:
659+
if initial_conformation.shape[0] != self.N:
660+
sys.stdout.write('Number of particles is not correct')
661+
sys.exit(0)
662+
if initial_conformation.shape[1] != 3:
663+
sys.stdout.write('The dimension should be three')
664+
sys.exit(0)
665+
self.xyz = initial_conformation
666+
667+
self.traj = []
668+
for t in tqdm(range(T)):
669+
if t % update == 0:
670+
self.updateXYZ()
671+
if t % every == 0:
672+
self.updateXYZ()
673+
self.traj.append(self.xyz)
674+
#sys.stdout.write('\rTimestep {}'.format(t+1))
675+
#sys.stdout.flush()
676+
self.updateModes(method=method)
677+
678+
self.traj = np.array(self.traj)
679+
680+
def reset(self):
681+
self.generateXYZ()
682+
448683

449684
@click.command()
450685
@click.argument('input', nargs=1)

0 commit comments

Comments
 (0)