-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathelastic.py
370 lines (290 loc) · 13.1 KB
/
elastic.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
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
import numpy as np
import numpy.linalg as la
from scipy.constants import elementary_charge
def convert_to_GPa(x):
""" convert eV/A^3 to GPa """
return x*elementary_charge*(10**10)**3/(10**9)
def convert_from_GPa(x):
""" convert GPa to eV/A^3 """
return x/(elementary_charge*(10**10)**3/(10**9))
def construct_C(crystalclass,Cijs):
"""
given the crystal class and elastic constants, construct the elastic stiffness tensor
I followed the same numbering system as in elastic.H in the anisotropic elasticity codes
** I have currently only implemented this for cubic (4), hexagonal (9) and isotropic (10) !!
Parameters
----------
crystalclass : 0 = triclinic (a != b != c, alpha!=beta!=gamma)
1 = monoclinic, diad || x_2 (a != b != c, alpha==gamma==90!=beta)
2: monoclinic, diad || x_3 (a != b != c, alpha==beta==90!=gamma)
3: orthorhombic (a != b != c, alpha==beta==gamma==90)
4: cubic _ (a == b == c, alpha==beta==gamma==90)
5: tetragonal (4 4 4|m) _ (a == b != c, alpha==beta==gamma==90)
6: tetragonal (4mm _422 4|mmm 42m) (a == b != c, alpha==beta==gamma==90)
7: trigonal, (3 3) _ (a == b != c, alpha==beta==90, gamma==120)
8: trigonal, (32 3m 3m) (a == b != c, alpha==beta==90, gamma==120)
9: hexagonal (a == b != c, alpha==beta==90, gamma==120)
10: isotropic
Cijs : list of elastic constants in GPa
Returns
-------
C : stiffness tensor C in Voigt notation, cartesian basis, GPa
"""
if crystalclass == 10: ## isotropic
C11,C12 = Cijs
C = np.array([[C11,C12,C12,0,0,0],
[C12,C11,C12,0,0,0],
[C12,C12,C11,0,0,0],
[0,0,0,(C11-C12)/2.,0,0],
[0,0,0,0,(C11-C12)/2.,0],
[0,0,0,0,0,(C11-C12)/2.]])
elif crystalclass == 4: ## cubic
C11,C12,C44 = Cijs
C = np.array([[C11,C12,C12,0,0,0],
[C12,C11,C12,0,0,0],
[C12,C12,C11,0,0,0],
[0,0,0,C44,0,0],
[0,0,0,0,C44,0],
[0,0,0,0,0,C44]])
elif crystalclass == 9: ## hexagonal
C11,C12,C13,C33,C44 = Cijs
C = np.array([[C11,C12,C13,0,0,0],
[C12,C11,C13,0,0,0],
[C13,C13,C33,0,0,0],
[0,0,0,C44,0,0],
[0,0,0,0,C44,0],
[0,0,0,0,0,(C11-C12)/2.]])
else:
raise ValueError('invalid crystal class!')
return C
def expand_C(C_voigt):
"""
expand the 6x6 C in Voigt notation to full 3x3x3x3 C tensor
Parameters
----------
C_voigt : stiffness tensor C in Voight notation, in cartesian basis
Returns
-------
C : 3x3x3x3 stiffness tensor C in cartesian basis
"""
## Voight notation: xx, yy, zz, yz, xz, xy
mapping = [[0,0],[1,1],[2,2],[1,2],[0,2],[0,1]]
C = np.zeros((3,3,3,3))
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
C[i,j,k,l] = C_voigt[mapping.index([min(i,j),max(i,j)]),
mapping.index([min(k,l),max(k,l)])]
return C
def rotate_C(C,M):
"""
given the stiffness tensor C in cubic basis (pqrs indices),
rotate it to mnt basis (ijkl indices)
C_ijkl = R_ip.R_jq.R_kr.R_ls.C_pqrs
Parameters
----------
C : 3x3x3x3 stiffness tensor C in cartesian basis
M : 3x3 matrix for rotating from mnt basis to cartesian basis
Returns
-------
C_mnt : 3x3x3x3 stiffness tensor C in mnt basis
"""
R = M.T
C_mnt = np.tensordot(R,np.tensordot(R,np.tensordot(R,np.tensordot(R,C,
axes=([1],[0])),axes=([1],[1])),axes=([1],[2])),axes=([1],[3]))
return C_mnt
#def Lambda2(k,RD,Dflat,M,A,a0):
#
# """
# Calculates Lambda2(k) = leading term in D(k)
# Lambda2(k) = -0.5 * summation_over_R(D(R)*(k_hat.R)^2)
#
# ** THIS VERSION BASED ON FORCE-CONSTANTS INSTEAD OF ELASTIC CONSTANTS
# IS NOT APPLICABLE TO MULTIPLE ATOM BASIS AND IS DEPRECIATED!! **
#
# Parameters
# ----------
# k : k_hat vector (direction only, magnitude 1)
# RD : list of vectors between pairs of atoms within the thin slab,
# given in primitive cell basis
# Dflat : list of "flattened" 3x3 force constant matrices corresponding to the atom pairs,
# listed in the same order as RD
# M : 3x3 matrix for rotating from mnt basis to cartesian basis
# A : 3x3 matrix for rotating from primitive cell basis to cartesian basis
# a0 : lattice constant in angstroms
#
# Returns
# -------
# L2 : 3x3 Lambda2(k) matrix evaluated for the given k vector
#
# Notes
# -----
# Refer to: D.R. Trinkle, Phys. Rev. B 78, 014110 (2008)
#
# """
#
# L2 = np.zeros((3,3))
#
# for RDi,Di in zip(RD,Dflat):
# ## convert RD vectors to mnt basis for this
# L2 -= 0.5 * Di * (np.dot(k,a0*np.dot(M.T,np.dot(A,RDi))))**2
#
# return L2
#def EGF_Fcoeffs(N,RD,Dflat,M,A,a0):
#
# """
# Calculates EGF(k) = inv(D(k)) = inv(Lambda2(k)),
# then expands EGF as a fourier series in the mn plane.
#
# ** THIS VERSION BASED ON FORCE-CONSTANTS INSTEAD OF ELASTIC CONSTANTS
# IS NOT APPLICABLE TO MULTIPLE ATOM BASIS AND IS DEPRECIATED!! **
#
# Parameters
# ----------
# N : number of k vectors to generate and compute EGF(k) for
# RD : list of vectors between pairs of atoms within the thin slab,
# given in primitive cell basis
# Dflat : list of "flattened" 3x3 force constant matrices corresponding to the atom pairs,
# listed in the same order as RD
# M : 3x3 matrix for rotating from mnt basis to cartesian basis
# A : 3x3 matrix for rotating from primitive cell basis to cartesian basis
# a0 : lattice constant in angstroms
#
# Returns
# -------
# GEn : list of complex ndarrays of fourier coefficients
# evaluated for each component of the EGF, i.e.
# [coeffs for xx component, coeffs for yy component, coeffs for zz component,
# coeffs for yz component, coeffs for xz component, coeffs for xy component]
#
# Notes
# -----
# Refer to: D.R. Trinkle, Phys. Rev. B 78, 014110 (2008)
#
# """
#
# ## generate list of k vectors that lie on m-n plane; angle is wrt +m axis
# ## calculates Lambda2(k) = leading term in D(k), then inverts Lambda2(k) to get EGF(k)
# EGF = [la.inv(Lambda2(k,RD,Dflat,M,A,a0))
# for k in [[np.cos(2*np.pi*n/N),np.sin(2*np.pi*n/N),0] for n in range(N)]]
#
# ## create separate lists for each independent component of the 3x3 EGF matrices
# EGF_xx,EGF_xy,EGF_xz = list(zip(*(list(zip(*EGF))[0])))
# EGF_yx,EGF_yy,EGF_yz = list(zip(*(list(zip(*EGF))[1])))
# EGF_zx,EGF_zy,EGF_zz = list(zip(*(list(zip(*EGF))[2])))
#
# return np.fft.fft([EGF_xx,EGF_yy,EGF_zz,EGF_yz,EGF_xz,EGF_xy])/len(EGF)
def Lambda2_kCk(k,C,M):
"""
compute Lambda2(k) = V*[k.C.k]
[k.C.k]_qr = sum_p,s (k[p]*C_pqrs*k[s])
Parameters
----------
k : unit vector in the mn plane, expressed in mnt coords
C : 3x3x3x3 stiffness tensor C in cartesian basis
M : 3x3 matrix for rotating from mnt basis to cartesian basis
Returns
-------
L2 : 3x3 Lambda2(k) matrix evaluated for the given k vector
"""
k = np.dot(M,k) ## convert k in the mn plane to cartesian coords
L2 = np.tensordot(k,np.tensordot(C,k,axes=([3],[0])),axes=([0],[0]))
## I've ignored the V here as it will cancel out later...
return L2
def EGF_Fcoeffs(N,C,M):
"""
Calculates EGF(k) = inv(D(k)) = inv(Lambda2(k)),
then expands EGF as a fourier series in the mn plane.
Parameters
----------
N : number of k vectors to generate and compute EGF(k) for
C : 3x3x3x3 stiffness tensor C in cartesian basis
M : 3x3 matrix for rotating from mnt basis to cartesian basis
Returns
-------
GEn : list of complex ndarrays of fourier coefficients
evaluated for each component of the EGF, i.e.
[coeffs for xx component, coeffs for yy component, coeffs for zz component,
coeffs for yz component, coeffs for xz component, coeffs for xy component]
Notes
-----
Refer to: D.R. Trinkle, Phys. Rev. B 78, 014110 (2008)
"""
## generate list of k vectors that lie on m-n plane; angle is wrt +m axis
## calculates Lambda2(k) = leading term in D(k), then inverts Lambda2(k) to get EGF(k)
EGF = [la.inv(Lambda2_kCk(k,C,M))
for k in [[np.cos(2*np.pi*n/N),np.sin(2*np.pi*n/N),0] for n in range(N)]]
## create separate lists for each independent component of the 3x3 EGF matrices
EGF_xx,EGF_xy,EGF_xz = list(zip(*(list(zip(*EGF))[0])))
EGF_yx,EGF_yy,EGF_yz = list(zip(*(list(zip(*EGF))[1])))
EGF_zx,EGF_zy,EGF_zz = list(zip(*(list(zip(*EGF))[2])))
return np.fft.fft([EGF_xx,EGF_yy,EGF_zz,EGF_yz,EGF_xz,EGF_xy])/len(EGF)
def G_largeR_ang(GEn,N,N_max):
"""
Calculates the angular term in the real space large R LGF
for a given finite number of angular values only.
Parameters
----------
GEn : list of complex ndarrays of fourier coefficients evaluated for each component of EGF
N : total number of fourier components in the fourier series;
also equal to the total number of angular values for which the
angular term in the real space large R LGF is to be explicitly computed
N_max : maximum number of fourier components to consider in the truncated fourier series
Returns
-------
phi_R_grid : list of angular terms in the real space large R LGF,
computed for N equally-spaced angular (phi) values.
Same as GEn, each entry in the list corresponds to
the values for the different components of the LGF:
[xx components, yy comps, zz comps, yz comps, xz comps, xy comps]
Notes
-----
Refer to: D.R. Trinkle, Phys. Rev. B 78, 014110 (2008)
"""
ang_coeffs = np.zeros([6,N],dtype=np.complex128)
for ang_coeffs_comp,GEn_comp in zip(ang_coeffs,GEn):
## evaluate separately for each of the 6 independent components
for n in range (1,N_max+1):
## include coeffs of exp(i*n*phi) as well as exp(-i*n*phi) in the truncated series
ang_coeffs_comp[n] += (((-1.)**(n/2))/n)*GEn_comp[n]
ang_coeffs_comp[-n] += (((-1.)**(n/2))/n)*GEn_comp[-n]
return np.fft.ifft(ang_coeffs)*N
def G_largeR(GEn,phi_R_grid,R,phi,N,t_mag):
"""
Calculates the real space large R LGF for given R, phi values.
Parameters
----------
GEn : list of complex ndarrays of fourier coefficients evaluated for each component of EGF
phi_R_grid : list of angular terms in the real space large R LGF, computed for
N equally-spaced angular (phi) values. Same as GEn, each entry in
the list corresponds to the values for the different components of LGF.
R : in-plane distance (on the m-n plane) between 2 atoms
phi : in-plane angle between 2 atoms, measured relative to +m axis
N : number of angular values for which the angular term in the real space large R LGF
has been explicitly computed
t_mag : magnitude of the periodic vector along the dislocation threading direction
(in Angstroms)
Returns
-------
G_largeR : 3x3 real space large R LGF matrix between 2 points
separated by distance R and angle phi
Notes
-----
Refer to: D.R. Trinkle, Phys. Rev. B 78, 014110 (2008)
"""
phi_ongrid = phi/(2*np.pi/N)
## angular values between which to interpolate
lower = int(np.floor(phi_ongrid))
upper = int(lower + 1)
lever_lower = phi_ongrid-lower
lever_upper = upper-phi_ongrid
logR = np.log(R)
## evaluate separately for each of the 6 independent components
## I've ignored the V here as it cancels out the V in the Lambda(2) term
G_largeR_comps = [((1./(2*np.pi*t_mag)) * (-GEn_comp[0]*logR +
lever_lower*phi_R_grid_comp[upper%N] + lever_upper*phi_R_grid_comp[lower]))
for GEn_comp,phi_R_grid_comp in zip(GEn,phi_R_grid)]
return [[G_largeR_comps[0].real,G_largeR_comps[5].real,G_largeR_comps[4].real],
[G_largeR_comps[5].real,G_largeR_comps[1].real,G_largeR_comps[3].real],
[G_largeR_comps[4].real,G_largeR_comps[3].real,G_largeR_comps[2].real]]