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