RMF
motion.py

Do a silly little simulation to show things moving around

1 ## \example motion.py
2 # Do a silly little simulation to show things moving around
3 
4 from __future__ import print_function
5 import RMF
6 import random
7 
8 number_of_particles = 10
9 number_of_frames = 100
10 minimum_radius = .5
11 maximum_radius = 3
12 box_size = 15
13 
14 # find the name for a temporary file to use to for writing the hdf5 file
15 tfn = RMF._get_temporary_file_path("motion.rmf3")
16 print("File is", tfn)
17 
18 # open the temporary file, clearing any existing contents
19 fh = RMF.create_rmf_file(tfn)
20 
21 tf = RMF.TypedFactory(fh)
22 pf = RMF.ParticleFactory(fh)
23 bf = RMF.BondFactory(fh)
24 sf = RMF.SegmentFactory(fh)
25 rf = RMF.RepresentationFactory(fh)
26 scf = RMF.ScoreFactory(fh)
27 
28 coords = []
29 radii = []
30 particles = []
31 
32 
33 def get_distance(v0, r0, v1, r1):
34  return (
35  ((v0[0] - v1[0]) ** 2 + (v0[1] - v1[1])
36  ** 2 + (v0[1] - v1[1]) ** 2) ** .5 - r0 - r1
37  )
38 
39 
40 def intersects(c, r, cs, rs):
41  for oc, orad in zip(cs, rs):
42  if get_distance(c, r, oc, orad) < 0:
43  # print c,r, "intersects", oc, orad
44  return True
45  return False
46 
47 
48 def too_far(c, r, oc, orad):
49  d = get_distance(c, r, oc, orad)
50  if d > .5:
51  return True
52  return False
53 
54 
55 def random_coordinates():
56  return RMF.Vector3(random.uniform(0, box_size),
57  random.uniform(0, box_size),
58  random.uniform(0, box_size))
59 
60 for p in range(0, number_of_particles):
61  cur = fh.get_root_node().add_child(str(p), RMF.REPRESENTATION)
62  tf.get(cur).set_type_name(str(p % 3))
63  pd = pf.get(cur)
64  pd.set_static_radius(random.uniform(minimum_radius, maximum_radius))
65  pd.set_static_mass(pd.get_radius())
66  curcoord = random_coordinates()
67  while intersects(curcoord, pd.get_radius(), coords, radii):
68  curcoord = random_coordinates()
69  radii.append(pd.get_radius())
70  coords.append(curcoord)
71  particles.append(cur)
72 
73 bn = fh.get_root_node().add_child("bonds", RMF.ORGANIZATIONAL)
74 
75 bonds = []
76 print("creating bonds")
77 for i in range(10):
78  bd0 = random.choice(particles)
79  bd1 = random.choice(particles)
80  if bd0 == bd1:
81  continue
82  bp = (particles.index(bd0), particles.index(bd1))
83  tries = 0
84  while too_far(coords[bp[0]], radii[bp[0]],
85  coords[bp[1]], radii[bp[1]])\
86  or intersects(coords[bp[0]], radii[bp[0]],
87  coords[:bp[0]] + coords[bp[0] + 1:],
88  radii[:bp[0]] + radii[bp[0] + 1:]):
89  coords[bp[0]] = random_coordinates()
90  tries = tries + 1
91  if tries == 100:
92  break
93  if tries < 100:
94  b = bn.add_child("bond", RMF.BOND)
95  bd = bf.get(b)
96  bd.set_static_bonded_0(bd0.get_id().get_index())
97  bd.set_static_bonded_1(bd1.get_id().get_index())
98 
99  bonds.append(bp)
100  coords[bonds[-1][0]] = coords[bonds[-1][1]]
101 
102 print("created", len(bonds), "bonds")
103 fn = fh.get_root_node().add_child("restraints", RMF.ORGANIZATIONAL)
104 features = []
105 feature_nodes = []
106 print("creating features")
107 for i in range(4):
108  f0 = random.choice(particles)
109  f1 = random.choice(particles)
110  if f0 == f1:
111  continue
112  f = fn.add_child("r", RMF.FEATURE)
113  rd = rf.get(f)
114  rd.set_representation([f0, f1])
115  features.append((particles.index(f0), particles.index(f1)))
116  feature_nodes.append(f)
117 
118 
119 def get_particle_distance(p0, p1):
120  pd0 = pf.get(p0)
121  pd1 = pf.get(p1)
122  v0 = pd0.get_coordinates()
123  v1 = pd1.get_coordinates()
124  r0 = pd0.get_radius()
125  r1 = pd1.get_radius()
126  return (
127  ((v0[0] - v1[0]) ** 2 + (v0[1] - v1[1])
128  ** 2 + (v0[1] - v1[1]) ** 2) ** .5 - r0 - r1
129  )
130 
131 bn = fh.get_root_node().add_child("box", RMF.ORGANIZATIONAL)
132 # draw box
133 bs = box_size
134 for e in [(RMF.Vector3(0, 0, 0), RMF.Vector3(0, 0, bs)),
135  (RMF.Vector3(0, 0, 0), RMF.Vector3(0, bs, 0)),
136  (RMF.Vector3(0, 0, 0), RMF.Vector3(bs, 0, 0)),
137  (RMF.Vector3(bs, bs, bs), RMF.Vector3(0, bs, bs)),
138  (RMF.Vector3(bs, bs, bs), RMF.Vector3(bs, 0, bs)),
139  (RMF.Vector3(bs, bs, bs), RMF.Vector3(bs, bs, 0)),
140  (RMF.Vector3(0, 0, bs), RMF.Vector3(0, bs, bs)),
141  (RMF.Vector3(0, 0, bs), RMF.Vector3(bs, 0, bs)),
142  (RMF.Vector3(0, bs, 0), RMF.Vector3(0, bs, bs)),
143  (RMF.Vector3(0, bs, 0), RMF.Vector3(bs, bs, 0)),
144  (RMF.Vector3(bs, 0, 0), RMF.Vector3(bs, bs, 0)),
145  (RMF.Vector3(bs, 0, 0), RMF.Vector3(bs, 0, bs))
146  ]:
147  en = bn.add_child("edge", RMF.GEOMETRY)
148  sd = sf.get(en)
149  sd.set_static_coordinates_list(e)
150 
151 
152 def write_frame(name):
153  fh.add_frame(name, RMF.FRAME)
154  for c, p in zip(coords, particles):
155  pd = pf.get(p)
156  pd.set_frame_coordinates(c)
157  for f, p in zip(feature_nodes, features):
158  s = get_particle_distance(particles[p[0]], particles[p[1]])
159  sd = scf.get(f)
160  sd.set_frame_score(s)
161 
162 write_frame("initial")
163 
164 for i in range(0, number_of_frames):
165  for i in range(len(coords)):
166  c = coords[i]
167  cn = RMF.Vector3(c[0] + random.uniform(-.2, .2),
168  c[1] + random.uniform(-.2, .2),
169  c[2] + random.uniform(-.2, .2))
170  if intersects(cn, radii[i], coords[:i] + coords[i + 1:], radii[:i] + radii[i + 1:]):
171  continue
172  bad = False
173  for b in [b for b in bonds if b[0] == i or b[1] == i]:
174  if b[1] == i:
175  b = [b[1], b[0]]
176  if too_far(cn, radii[i], coords[b[1]], radii[b[1]]):
177  bad = False
178  break
179  coords[i] = cn
180  write_frame(str(i))