# EventDrivenMD / edmd.py

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237``` ```import numpy as np import math import itertools class Particle(object): def __init__(self, x, v): self.x = np.array(x) self.v = np.array(v) self.event_times = [] def __repr__(self): return "Particle(%s, %s)" % (self.x, self.v) def move(self, delta_t): self.x += delta_t * self.v def overlaps_with(self, other_particle, radius): squared_distance = np.sum((other_particle.x - self.x) ** 2) return squared_distance < 4.0 * radius ** 2 class Box(object): def __init__(self, size=10): self.size = size def contains(self, position, radius): return position[0] <= self.size - radius and \ position[0] >= -self.size + radius and \ position[1] <= self.size - radius and \ position[1] >= -self.size + radius class CollisionRules(object): def __init__(self, box=None, eps=1.0, radius=1.0): self.tolerance = 0.0e-10 self.eps = eps self.radius = radius self.box = box if box else Box() def collide(self, particle_i, particle_j): x_ij = particle_i.x - particle_j.x v_ij = particle_i.v - particle_j.v unit_vector = x_ij / np.linalg.norm(x_ij) h = (1 + self.eps) * np.dot(v_ij, unit_vector) / 2.0 particle_i.v -= h * unit_vector particle_j.v += h * unit_vector def collide_with_wall(self, particle): if particle.x[1] < 0: if particle.x[1] > particle.x[0] or particle.x[1] > -particle.x[0]: particle.v[0] = -particle.v[0] else: particle.v[1] = -particle.v[1] else: if particle.x[1] < -particle.x[0] or particle.x[1] < particle.x[0]: particle.v[0] = -particle.v[0] else: particle.v[1] = -particle.v[1] def time_to_next_collision(self, particle_i, particle_j): x_ij = particle_i.x - particle_j.x v_ij = particle_i.v - particle_j.v x_ij_dot_v_ij = np.dot(x_ij, v_ij) if x_ij_dot_v_ij > 0: return float('inf') squared_distance = np.sum(x_ij ** 2.0) # all particles with same radius if squared_distance >= 4.0 * (self.radius ** 2): radius = self.radius else: radius = self.radius - self.tolerance q = squared_distance - 4.0 * radius ** 2 squared_velociy_difference = np.sum(v_ij ** 2) w = x_ij_dot_v_ij ** 2 - q * squared_velociy_difference if w < 0: return float('inf') delta_t = q / (-x_ij_dot_v_ij + math.sqrt(w)) x_i = particle_i.x + delta_t * particle_i.v if not self.box.contains(x_i, self.radius): return float('inf') x_j = particle_j.x + delta_t * particle_j.v if not self.box.contains(x_j, self.radius): return float('inf') return delta_t def time_to_next_collision_with_wall(self, part): radius = self.radius if any(np.abs(part.x - self.box.size + self.radius) > 0): radius -= self.tolerance result = min(self._compute_dt(part, 0), self._compute_dt(part, 1)) return result def _compute_dt(self, particle, index): dtime = float('inf') if particle.v[index] > 0: dtime = (self.box.size - (particle.x[index] + self.radius)) /\ particle.v[index] elif particle.v[index] < 0: dtime = (self.box.size + (particle.x[index] - self.radius)) /\ -particle.v[index] return dtime class Simulation(object): def __init__(self, particles, collision_rules=None, nprint=1): self.particles = particles if collision_rules is None: collision_rules = CollisionRules() self.collision_rules = collision_rules self.events = dict() self.time = 0 self.update_particles() self.nprint = nprint self.t_events = [] def update_particles(self): for particle in self.particles: self.update(particle) def kinenergy(self): return sum(np.sum(particle.v ** 2) for particle in self.particles) def update(self, particle): for time in particle.event_times: if time in self.events: del self.events[time] particle.event_times = [] for other_particle in self.particles: if other_particle != particle: dt = self.collision_rules.time_to_next_collision(particle, other_particle) if dt < float('inf'): time_coll = self.time + dt self.events[time_coll] = (particle, other_particle) particle.event_times.append(time_coll) other_particle.event_times.append(time_coll) dt = self.collision_rules.time_to_next_collision_with_wall(particle) # if dt >= 0: if dt < float('inf'): time_coll = self.time + dt self.events[time_coll] = (particle, None) particle.event_times.append(time_coll) def move_particles(self, delta_t): for particle in self.particles: particle.move(delta_t) def run(self, steps=10): for i in xrange(steps): t_next_event = min(self.events.keys()) next_event = self.events[t_next_event] if i % self.nprint == 0: print i, self.time self.t_events.append(self.time) if next_event[1] is None: delta_t = self.collision_rules.time_to_next_collision_with_wall(next_event[0]) self.move_particles(delta_t) self.time = t_next_event self.collision_rules.collide_with_wall(next_event[0]) self.update(next_event[0]) else: delta_t = self.collision_rules.time_to_next_collision(next_event[0], next_event[1]) self.move_particles(delta_t) self.time = t_next_event self.collision_rules.collide(next_event[0], next_event[1]) self.update(next_event[0]) self.update(next_event[1]) def load_initial_particles(coordinates_file_path, velocities_file_path): coords = np.loadtxt(coordinates_file_path) vels = np.loadtxt(velocities_file_path) particles = [Particle(x, v) for x, v in itertools.izip(coords, vels)] return particles def create_random_particles(seed, box, num_particles, radius, minimium_velocity, maximum_velocity): particles = [] np.random.seed(seed) coords_upper_bound = box.size - radius coords_lower_bound = -box.size + radius for i in xrange(0, num_particles): overlap = True while(overlap): x = np.random.uniform(coords_lower_bound, coords_upper_bound, 2) v = np.random.uniform(minimium_velocity, maximum_velocity, 2) new_particle = Particle(x, v) overlap = is_there_overlap(new_particle, particles, radius) particles.append(new_particle) return particles def is_there_overlap(new_particle, particles, radius): for particle in particles: if new_particle.overlaps_with(particle, radius): return True return False def create_simulation_from_files(coords_filename, vels_filename, box_size, eps): particles = load_initial_particles(coords_filename, vels_filename) rules = CollisionRules(Box(box_size), eps) return Simulation(particles, rules) def create_simulation(num_particles, seed, radius, box_size, eps, minimium_velocity, maximum_velocity): box = Box(box_size) particles = create_random_particles(seed, box, num_particles, radius, minimium_velocity, maximum_velocity) rules = CollisionRules(box, eps, radius) return Simulation(particles, rules) if __name__ == '__main__': import sys coords_filename = "./tests/data/coords_ini_test_update.txt" if len(sys.argv) > 1: coords_filename = sys.argv[1] vels_filename = "./tests/data/v_ini_test_update.txt" if len(sys.argv) > 2: vels_filename = sys.argv[2] sim = create_simulation_from_files(coords_filename, vels_filename, box_size=10) sim.run(100) ```