Coverage for gamdpy/simulation_boxes/lees_edwards.py: 87%
142 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4Created on Thu Jan 30 13:49:08 2025
6@author: nbailey
7"""
9import numpy as np
10import numba
11from numba import cuda
12import math
13from .simulationbox import SimulationBox
16class LeesEdwards(SimulationBox):
17 """ Simulation box class with LeesEdwards bondary conditions
19 Example
20 -------
22 >>> import gamdpy as gp
23 >>> import numpy as np
24 >>> simbox = gp.LeesEdwards(D=3, lengths=np.array([3,4,5]), box_shift=1.0)
26 """
27 def __init__(self, D, lengths, box_shift=0., box_shift_image=0):
28 if D < 2:
29 raise ValueError("Cannot use LeesEdwards with dimension smaller than 2")
30 self.D = D
31 self.box_shift = box_shift
32 self.box_shift_image = np.float32(box_shift_image)
33 self.data_array = np.zeros(D+2, dtype=np.float32)
34 self.data_array[:D] = np.array(lengths, dtype=np.float32)
35 self.data_array[D] = self.box_shift
36 self.data_array[D+1] = self.box_shift_image
37 self.len_sim_box_data = D+2
39 return
41 def get_name(self):
42 return "LeesEdwards"
44 def copy_to_device(self):
45 self.d_data = cuda.to_device(self.data_array)
48 def copy_to_host(self):
49 D = self.D
50 self.data_array = self.d_data.copy_to_host()
51 self.box_shift = self.data_array[D]
52 self.boxshift_image = self.data_array[D+1]
54 def get_volume_function(self):
55 D = self.D
56 def volume(sim_box):
57 ''' Returns volume of the rectangular box '''
58 vol = sim_box[0]
59 for i in range(1,D):
60 vol *= sim_box[i]
61 return vol
62 return volume
64 def get_lengths(self):
65 return self.data_array[:self.D].copy()
67 def get_volume(self):
68 #self.copy_to_host()
69 return self.get_volume_function()(self.data_array[:self.D])
71 def scale(self, scale_factor):
72 self.data_array[:self.D] *= scale_factor
75 def get_dist_sq_dr_function(self):
76 """Generates function dist_sq_dr which computes displacement and distance for one neighbor """
77 D = self.D
79 def dist_sq_dr_function(ri, rj, sim_box, dr):
80 ''' Returns the squared distance between ri and rj applying MIC and saves ri-rj in dr '''
81 box_shift = sim_box[D]
82 for k in range(D):
83 dr[k] = ri[k] - rj[k]
85 dist_sq = numba.float32(0.0)
86 box_1 = sim_box[1]
87 dr[0] += (-box_shift if numba.float32(2.0) * dr[1] > +box_1 else
88 (+box_shift if numba.float32(2.0) * dr[1] < -box_1 else
89 numba.float32(0.0)))
91 for k in range(D):
92 box_k = sim_box[k]
93 dr[k] += (-box_k if numba.float32(2.0) * dr[k] > +box_k else
94 (+box_k if numba.float32(2.0) * dr[k] < -box_k else numba.float32(0.0))) # MIC
95 dist_sq = dist_sq + dr[k] * dr[k]
96 return dist_sq
98 return dist_sq_dr_function
100 def get_dist_sq_function(self):
102 D = self.D
103 def dist_sq_function(ri, rj, sim_box):
104 ''' Returns the squared distance between ri and rj applying MIC'''
105 box_shift = sim_box[D]
106 dist_sq = numba.float32(0.0)
108 # first shift the x-component depending on whether the y-component is wrapped
109 dr1 = ri[1] - rj[1]
110 box_1 = sim_box[1]
111 x_shift = (-box_shift if numba.float32(2.0) * dr1 > box_1 else
112 (+box_shift if numba.float32(2.0) * dr1 < -box_1 else
113 numba.float32(0.0)))
114 # then wrap as usual for all components
115 for k in range(D):
116 dr_k = ri[k] - rj[k]
117 if k == 0:
118 dr_k += x_shift
119 box_k = sim_box[k]
120 dr_k += (-box_k if numba.float32(2.0) * dr_k > +box_k else
121 (+box_k if numba.float32(2.0) * dr_k < -box_k else numba.float32(0.0))) # MIC
122 dist_sq = dist_sq + dr_k * dr_k
123 return dist_sq
124 return dist_sq_function
126 def get_apply_PBC(self):
127 D = self.D
128 def apply_PBC(r, image, sim_box):
130 # first shift the x-component depending on whether the y-component is outside the box
131 # note: assumes at most one box length needs to be added/subtracted.
132 box_shift, bs_image = sim_box[D], int(sim_box[D+1])
133 box1_half = sim_box[1] * numba.float32(0.5)
134 if r[1] > + box1_half:
135 r[0] -= box_shift
136 image[0] -= bs_image
137 if r[1] < -box1_half:
138 r[0] += box_shift
139 image[0] += bs_image
140 # then put everything back in the box as usual
141 for k in range(D):
142 if r[k] * numba.float32(2.0) > +sim_box[k]:
143 r[k] -= sim_box[k]
144 image[k] += 1
145 if r[k] * numba.float32(2.0) < -sim_box[k]:
146 r[k] += sim_box[k]
147 image[k] -= 1
148 return
149 return apply_PBC
151 def get_update_box_shift(self):
152 D = self.D
153 def update_box_shift(sim_box, shift):
154 # carry out the addition in double precision
155 sim_box[D] = numba.float32(sim_box[D] + numba.float64(shift))
156 Lx = sim_box[0]
157 Lx_half = Lx*numba.float32(0.5)
158 if sim_box[D] > +Lx_half:
159 sim_box[D] -= Lx
160 sim_box[D+1] += 1
161 if sim_box[D] < -Lx_half:
162 sim_box[D] += Lx
163 sim_box[D+1] -= 1
164 return
165 return update_box_shift
167 #def get_dist_moved_sq_function(self):
168 # D = self.D
169 # def dist_moved_sq_function(r_current, r_last, sim_box, sim_box_last):
170 # zero = numba.float32(0.)
171 # half = numba.float32(0.5)
172 # one = numba.float32(1.0)
173 # box_shift = sim_box[D]
174 # dist_moved_sq = zero
177 # strain_change = sim_box[D] - sim_box_last[D] # change in box-shift
178 # strain_change += (sim_box[D+1] - sim_box_last[D+1]) * sim_box[0] # add contribution from box_shift_image
179 # strain_change /= sim_box[1] # convert to (xy) strain
182 # # we will shift the x-component when the y-component is 'wrapped'
183 # dr1 = r_current[1] - r_last[1]
184 # box_1 = sim_box[1]
185 # y_wrap = (one if dr1 > half*box_1 else
186 # -one if dr1 < -half*box_1 else zero)
188 # x_shift = y_wrap * box_shift + (r_current[1] -
189 # y_wrap*box_1) * strain_change
190 # # see the expression in Chatoraj Ph.D. thesis. Adjusted here to
191 # # take into account BC wrapping (otherwise would use the images
192 # # ie unwrapped positions)
194 # for k in range(D):
195 # dr_k = r_current[k] - r_last[k]
196 # if k == 0:
197 # dr_k -= x_shift
198 # box_k = sim_box[k]
199 # dr_k += (-box_k if numba.float32(2.0) * dr_k > +box_k else
200 # (+box_k if numba.float32(2.0) * dr_k < -box_k else numba.float32(0.0)))
201 # dist_moved_sq = dist_moved_sq + dr_k * dr_k
204 # return dist_moved_sq
205 # return dist_moved_sq_function
207 def get_dist_moved_exceeds_limit_function(self):
208 D = self.D
210 def dist_moved_exceeds_limit_function(r_current, r_last, sim_box, sim_box_last, skin, cut):
211 """
212 Returns true if the distance moved since last neighbor-list update exceeds half the skin, taking the change in box-shift into account
213 See Chattoraj PhD thesis for criterion for neighbor list checking under shear https://pastel.hal.science/pastel-00664392/
214 """
215 zero = numba.float32(0.)
216 half = numba.float32(0.5)
217 one = numba.float32(1.0)
218 box_shift = sim_box[D]
219 dist_moved_sq = zero
222 strain_change = sim_box[D] - sim_box_last[D] # change in box-shift
223 strain_change += (sim_box[D+1] - sim_box_last[D+1]) * sim_box[0] # add contribution from box_shift_image
224 strain_change /= sim_box[1] # convert to (xy) strain
226 # we will shift the x-component when the y-component is 'wrapped'
227 dr1 = r_current[1] - r_last[1]
228 box_1 = sim_box[1]
229 y_wrap = (one if dr1 > half*box_1 else
230 -one if dr1 < -half*box_1 else zero)
232 x_shift = y_wrap * box_shift + (r_current[1] -
233 y_wrap*box_1) * strain_change
234 # see the expression in Chatoraj Ph.D. thesis. Adjusted here to
235 # take into account BC wrapping (otherwise would use the images
236 # ie unwrapped positions)
238 for k in range(D):
239 dr_k = r_current[k] - r_last[k]
240 if k == 0:
241 dr_k -= x_shift
242 box_k = sim_box[k]
243 dr_k += (-box_k if numba.float32(2.0) * dr_k > +box_k else
244 (+box_k if numba.float32(2.0) * dr_k < -box_k else numba.float32(0.0)))
245 dist_moved_sq = dist_moved_sq + dr_k * dr_k
247 skin_corrected = skin - abs(strain_change)*cut
248 if skin_corrected < zero:
249 skin_corrected = zero
251 return dist_moved_sq > skin_corrected*skin_corrected*numba.float32(0.25)
253 return dist_moved_exceeds_limit_function
255 def get_loop_x_addition(self):
256 return 1
258 def get_loop_x_shift_function(self):
259 D = self.D
260 def loop_x_shift_function(sim_box, cell_length_x):
261 box_shift = sim_box[D]
262 return -int(math.ceil(box_shift/cell_length_x))
264 return loop_x_shift_function