-
Notifications
You must be signed in to change notification settings - Fork 4
/
Bravais_Lattice.m
389 lines (385 loc) · 19.7 KB
/
Bravais_Lattice.m
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
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
classdef Bravais_Lattice < handle
% Every operation is in accordance with the International tables of
% crystallography Vol A, A1
properties (SetAccess = public)
phase; % string specifying wheter parent/austenite or product/martensite
Bravais_type; % a string specifying the Bravais lattice type e.g. 'cubic'
Centering; % a string specifying the exact type e.g. simple 'P',
% face-centered 'F', body-centered 'I' etc.
% Lattice_parameters
Lp = [1.0, 1.0, 1.0, pi/2, pi/2, pi/2];
%
CPPs; % close packed planes of lattice
CP_dirs; % close packed directions of lattice
%
slip_plane_families; % this and the next variable are for writing the results
slip_dir_families;
slip_planes; % full array of slip planes created from plane_families
slip_directions; % full array of slip directions created from slip_dir_families
% the order of slip_planes and slip_direction fits together such
% that all independent deformations are taken care of (i.e. both deformation directions)
mirror_planes; % generally low_index planes
end
properties (Dependent, Access = public)
Point_group; % 3x3xN matrix containing the point group symmetry matrices
% This group is often also denoted as Laue-Group, which however
% also contains matrices generated in a left-handed system (det<0)
% which we do not consider here
E;
% primitive (crystallographic) basis vectors e_i stored as E(:,:,1) =
% [e1;e2;e3], row vectors
% i.e. lattice vectors with integer coefficients t_i: {t_vec = t1*a1 + t2*a2 + t3*a3}.
% this means that the a_i are translation vectors!
% This matrix array also contains the conversion matrix to the
% conventional basis in E(:,:,2)
C;
% conventional crystallographic = covariant = real basis = direct lattice
% c_i stored as C = [c1;c2;c3], (such that as many angles as
% possible amout to 90 degree. e.g.
% primitive lattice --> conventional lattice = primitive (e.g. primitive cubic),
% otherwise "centered" e.g. face centered cubic.
% note that for e.g. monoclinic, one conventional basis vector is
% not parrallel to a primitive basis vector
% These are the crystallographic bases used in the International Tables A !
% must be user specified because of different possibilities
% e.g. two for monoclinic
Cov_metric_C;
Contra_metric_C;
reciprocal_base; % = dual- or contravariant base
Lattice_group;
density;
%
end
methods
% constructor
function obj = Bravais_Lattice()
end
%------------------------------------------------------------------
function obj = set.Bravais_type(obj, input)
Bravais_strings = ['cubic '; 'hexagonal '; 'rhombohedral'; ...
'tetragonal '; 'tetragonal '; 'orthorhombic'; 'monoclinic '; ...
'triclinic '];
Bravais_list = cellstr(Bravais_strings);
if ismember( input, Bravais_list)
obj.Bravais_type = input;
else
error('Invalid Bravais Type')
end
end
%------------------------------------------------------------------
function obj = set.Centering(obj, input )
Centering_strings = ['P'; 'C'; 'I'; 'F'];
if ismember( input, Centering_strings)
obj.Centering = input;
else
error('Invalid Centering')
end
end
%------------------------------------------------------------------
% TODO write function or class for user input e.g.
% bravais_type = input('choose bravais lattice') %later maybe make a
% GUI for that#
function obj = set.Lp(obj, Lpars)
switch obj.Bravais_type % the number of lattice parameters must match the bravais type
case 'cubic'
obj.Lp(1:3) = Lpars;
% Par_num = [1, 0]; %number of different lengths and angles unequal 90�
% case 'hexagonal'
% Par_num = [2, 0];
% case 'monoclinic'
% Par_num = [3, 1];
% case 'triclinic'
% Par_num = [3, 3];
% %TODO implement other lattices if needed
end
% for i = 1 : obj.Par_num(2)
% val = input( ['choose angle(s) because the', ...
% 'assignment of an coordinate system is not unique: alpha:', ...
% 'angle between a-b, beta b-c, gamma c-a ']);
% if isnumeric(Lattice_params{i})
% obj.Lp{i} = val;
% else
% error('Value must be numeric')
% end
% deg2rad
% end
end
%------------------------------------------------------------------
function point_group = get.Point_group( obj )
point_group = matrix_group( obj.Bravais_type );
end
%------------------------------------------------------------------
function bool = in_PointGroup(g)
% Given another MatrixGroup g (3x3xN array), checks if its a subgoup of the
% lattice point group
% compare group order, i.e. number of elements in group
if ~ismember( size(g,3), divisors( size(Point_group,3)) )
bool = false;
return
else % check for matrices if they are contained
for i=1:size(g,3)
if ~in_matrix_array( g(:,:,3), Point_group )
bool = false;
return
end
end
bool = true;
end
end
%------------------------------------------------------------------
function value = get.E(obj)
% generate the base matrix "E" for primitive cell of 3D Bravais lattices
% it is stored in E(:,:,1)
% in E(:,:,2) the conversion matrix W from the primitive P to the
% conventional base C is stored like: C = E*W
% this is different to the internaltional tables Cryst. Vol A, p.81-83
% I took it from the python code pytrans from github and verified
% it for bcc to fcc.
switch obj.Bravais_type
case 'cubic'
switch obj.Centering
case 'P' % simple
value = obj.Lp(1) * eye(3);
case 'F' % Face centered fcc
value = 0.5 * obj.Lp(1)*[1, 1, 0;
0, 1, 1;
1, 0, 1]';
value(:,:,2) = 0.5* [1, 0, 1;
1, 1, 0;
0, 1, 1];
case 'I' % Body centered bcc
value = 0.5 * obj.Lp(1) * [1, 1, 1;
-1, 1, 1;
-1, -1, 1]';
value(:,:,2) = 0.5 * [1, -1, -1;
1, 1, -1;
1, 1, 1];
end
case 'hexagonal'
value = cat2(obj.Lp(1) *[1, 0, 0]', ...
obj.Lp(1)* [0.5, sqrt(3)/2, 0]', ...
0, obj.Lp(2)* [0, 0, 1]);
case 'trigonal' %<111> is the 3fold axis
% value =
c = cos(obj.Lp(1));
ty = sqrt((1. - c) / 6.);
tz = sqrt((1. + 2. * c) / 3.);
u = tz - 2. * sqrt(2.) * ty;
v = tz + sqrt(2.) * ty;
d = obj.Lp(1) / sqrt(3);
value = d*[u, v, v;
v, u, v;
v, v, u]';
case 'tetragonal'
if strcmpi(obj.Centering, 'P')
value = cat2(obj.Lp(1)* [1, 0, 0], ...
obj.Lp(1)* [0, 1, 0], ...
obj.Lp(3)* [0, 0, 1]);
elseif strcmpi(obj.Centering, 'I')
a = obj.Lp(1) / 2.;
value = [a, a, obj.Lp(3)/2.;
-a, a, obj.Lp(3)/2.;
-a, -a, obj.Lp(3)/2.]';
value(:,:,2) = 0.5 * [1, -1, -1;
1, 1, -1;
1, 1, 1];
end
case 'orthorhombic'
value = cat(2, obj.Lp(1)*[1 0 0]' , obj.Lp(2)*[0 1 0]' , obj.Lp(3)*[0 0 1]' );
% TODO
switch obj.Centering
case 'A' % 'B','C' one-face centred, % TODO
% # base centered orthorhombic
% a = p[0]
% b = p[1]
% c = p[2]
% e1 = np.array([a / 2, b / 2, 0])
% e2 = np.array([-a / 2, b / 2, 0])
% e3 = np.array([0, 0, c])
% C = np.array([[0.5, -0.5, 0],
% [0.5, 0.5, 0],
% [0, 0, 1]])
case 'F'
% # face centered orthorhombic
% a = p[0]
% b = p[1]
% c = p[2]
% e1 = np.array([a / 2, b / 2, 0])
% e2 = np.array([0, b / 2, c / 2])
% e3 = np.array([a / 2, 0, c / 2])
% C = 0.5 * np.array([[1, 0, 1],
% [1, 1, 0],
% [0, 1, 1]])
case 'I'
% # body centered orthorhombic
% a = p[0]
% b = p[1]
% c = p[2]
% e1 = np.array([a / 2, b / 2, c / 2])
% e2 = np.array([-a / 2, b / 2, c / 2])
% e3 = np.array([-a / 2, -b / 2, c / 2])
% C = 0.5 * np.array([[1, -1, -1],
% [1, 1, -1],
% [1, 1, 1]])
end
% add other monoclinic settings
case 'monoclinic' % unique axis b
switch obj.Centering
case 'P'
a = obj.Lp(1);
b = obj.Lp(2);
c = obj.Lp(3);
beta = obj.Lp(3);
value = cat(2, [a, 0, 0]', [0, b, 0]', [c*cos(beta), 0, c*sin(beta)]');
case 'A' % 'B','C' one-face centred, % TODO # base centered monoclinic
% a = p[0]
% b = p[1]
% c = p[2]
% beta = radians(p[3])
% e1 = np.array([a / 2, b / 2, 0])
% e2 = np.array([-a / 2, b / 2, 0])
% e3 = np.array([c * np.cos(beta), 0, c * np.sin(beta)])
% C = np.array([[0.5, -0.5, 0],
% [0.5, 0.5, 0],
% [0, 0, 1]])
end
% case 'triclinic'
% a = obj.Lp(1);
% b = obj.Lp(2);
% c = obj.Lp(3);
% alpha = radians(p[3])
% beta = radians(p[4])
% gamma = radians(p[5])
% e1 = np.array([a, 0, 0])
% e2 = np.array([b * np.cos(gamma), b * np.sin(gamma), 0])
% e3 = np.array([c * np.cos(beta), c * (np.cos(alpha) - np.cos(beta) * np.cos(gamma)) / np.sin(gamma),
% c * np.sqrt(1 + 2 * np.cos(alpha) * np.cos(beta) * np.cos(gamma) - np.cos(alpha) ** 2
% - np.cos(beta) ** 2 - np.cos(gamma) ** 2) / np.sin(gamma)])
%
% return np.array([e1, e2, e3]).T, la.inv(C)
end
if ~exist('value', 'var')
error('Bravais or centering type has not been set properly')
end
if size(value,3) == 1
value(:,:,2) = eye(3);
end
end
%------------------------------------------------------------------
function base = get.C(obj)
if size(obj.E,3) == 2
base = obj.E(:,:,1)* inverse( obj.E(:,:,2) );
else
base = obj.E(:,:,1);
end
end
%------------------------------------------------------------------
function metric = get.Cov_metric_C(obj)
% given the basis vectors E of the primitive crystallographic lattice
% calculate the corresponding covariant metric Z_i
metric = calc_metric(obj.E);
end
%------------------------------------------------------------------
function cont_metric = get.Contra_metric_C( obj )
% given the basis vectors E of the primitive crystallographic lattice
% calculate the corresponding Contravariant metric Z^i
% mf = matrix_functions;
cont_metric = inv(obj.Cov_metric_C);
end
%------------------------------------------------------------------
function db = get.reciprocal_base(obj)
% given the basis vectors E of the primitive crystallographic lattice
% calculate its dual basis
db = calc_dual_basis( obj.E );
end
%-------------------------------------------------------------------
% The components of any vector referred, or the coordinates of any point
% to the reciprocal basis represent the
% Miller indices of a plane whose normal is along that vector, with the spacing
% of the plane given by the inverse of the magnitude of that vector.
%%%%%
% at the moment the primitive lattice is of no concern
function prim = to_primitive_Miller( idx )
% convert Miller Indices or a list of Miller indices stored in
% an 3xN matrix from from conventional to primitive base
% base to a primitive base
for i = 1:size(idx,2)
prim(:,i) = C * idx(:,i);
end
end
%-------------------------------------------------------------------
function miller = to_conventional_Miller( idx )
% convert Miller indices or a list of Miller indices "idx"
% stored in a 3xN matrix from primitive to conventional base
% Miller indices are always integers!
%mf = matrix_functions();
for i = 1:size(idx,2)
miller(:,i) = inv( obj.E(:,:,2) )* idx(:,i);
end
end
%------------------------------------------------------------------
function obj = set.density( obj, Mz)
% Mz... Matrix holding in the first row M_i
% and in the second row z_i
% M...atomic Mole Masses, z...number of atom species in the unit cell
volume = sqrt(det(Cov_metric));
obj.density = (Mz(1,:).*Mz(2,:)) / volume * Avogadro;
end
%------------------------------------------------------------------
function vector = vec_from_coords(obj, coords )
% assemble a lattice vector based on its components or a point
% based on its coords
vector = coords(1)*obj.C(:,1) + coords(2)*obj.C(:,2) + coords(3)*obj.C(:,3);
end
%------------------------------------------------------------------
function length = vec_length( components )
% given the components of a lattice vector calculates its length
length = sqrt( components' * Cov_metric * components);
end
%------------------------------------------------------------------
function distance = distance_atoms(coords1, coords2)
X12 = coords1 - coords2;
distance = vec_length( X12 );
end
%------------------------------------------------------------------
function angle = bond_angle( coords1, coords2, coords3 )
% given three atom positions. calculates the atomic angle so
% thats it is measured at position of atom 1
X12 = coords1 - coords2;
X13 = coords1 - coords2;
a = X12' * Cov_metric * X13; % Z_ij = U^i V^i
angle = acos( a / (sqrt(a)* sqrt(X13' * Cov_metric * X12)));
end
%------------------------------------------------------------------
function components = components_tensorcalculus( vec, basis)
% given a vector and a co- or contravariant basis, this function yields the
% co- or contravariant components of the vector. (co -> co, contra-> contra)
% see Grinfeld - Intro to Tensorcalculus... p.61
for i = 1:size(basis,2)
components(i,1) = dot(vec, basis(:,i));
end
end
%------------------------------------------------------------------
% alternative solution of the above problem solving a coupled
% system of equations...
% function components = decomposition_wrt_to_basis( vec, basis )
% % Given a vector (colon) and a basis (matrix with colon vectors)
% % this function calculates the components of the vector w.r.t to the basis
% % (for affine bases only, i.e. constant basis in space)
% % See e.g. Grinfeld - Introduction to Tensor analysis calculus of moving surfaces - p61
% components = zeros(3,1);
% if isorthogonal(basis) % note: even simpler for an orthonormal basis
% for i=1:3
% components(i,1) = dot(vec , basis(:,i) ) / dot(basis(:,i),basis(:,i));
% end
% else % for affine bases (i.e constant in space but non-orthogonal)
% vec_dot_e_i = zeros(3,1);
% for i=1:3
% vec_dot_e_i(i,1) = dot(vec , basis(:,i) );
% end
% metric = metric_from_basis( basis );
% components(i,1) = metric \ vec_dot_e_i;
% end
% end
end % end methods
end % end class