wmaee.scopes.cij

  1"""
  2Convenient functions for working with elasticity.
  3"""
  4
  5import numpy as np
  6from typing import *
  7from ase import Atoms
  8from typing import Tuple
  9from typing import Union
 10from itertools import product
 11# from frozendict import frozendict as fdict
 12from numpy.typing import NDArray, ArrayLike
 13
 14from numpy import ndarray
 15
 16
 17
 18
 19
 20def get_ULICS(max_eps: float = 1.5e-2) -> ndarray:
 21    """
 22    Returns ULICS (universal linear-independent coupling strains) used
 23    for deformations in the stress-strain method.
 24
 25    For details, refer to doi:10.1016/j.cpc.2009.11.017.
 26
 27    Parameters
 28    ----------
 29    max_eps : float, optional
 30        Magnitude of the stress (largest component in Voigt's notation).
 31        Defaults to 1.5e-2.
 32
 33    Returns
 34    -------
 35    numpy.ndarray
 36        6x6 matrix of ULICS.
 37    """
 38    # Creating the 6x6 matrix of ULICS
 39    ULICS = max_eps / 6.0 * np.array([
 40        [1, -2, 3, -4, 5, -6],
 41        [2, 1, -5, -6, 4, 3],
 42        [3, 4, -1, 5, 6, -2],
 43        [4, -3, 6, 1, -2, 5],
 44        [5, 6, 2, -3, -1, -4],
 45        [6, -5, -4, 2, -3, 1]
 46    ], dtype=np.float64)
 47
 48    return ULICS
 49
 50
 51def apply_strain(struct: Atoms, strain: ArrayLike, div_two: bool = True) -> Atoms:
 52    """
 53    Applies strain to a structure.
 54
 55    Parameters
 56    ----------
 57    struct : ase.Atoms
 58        Initial structure.
 59    strain : ArrayLike
 60        Strain to be applied.
 61    div_two : bool, optional
 62        Whether to apply a factor of 2 when converting from
 63        Voigt's notation to tensorial 3x3 strain (default is True).
 64
 65    Returns
 66    -------
 67    ase.Atoms
 68        Deformed structure.
 69    """
 70    strain = np.array(strain)
 71
 72    # Try to convert the vector
 73    if strain.shape != (3, 3):
 74        strain = from_voigt(strain, div_two=div_two)
 75
 76    deformation_matrix = np.eye(3) + strain
 77
 78    new_cell = np.dot(deformation_matrix, struct.cell.array)
 79    new_struct = struct.copy()
 80    new_struct.set_cell(new_cell, scale_atoms=True)
 81
 82    return new_struct
 83
 84
 85
 86def index_from_voigt(n: int) -> Union[tuple, None]:
 87    """
 88    Convert Voigt's index to matrix indices.
 89
 90    Parameters
 91    ----------
 92    n : int
 93        Voigt's index.
 94
 95    Returns
 96    -------
 97    Union[tuple, None]
 98        Tuple representing matrix indices (a, b) if n is a valid Voigt's index, None otherwise.
 99    """
100    if n == 0:
101        return 0, 0
102    elif n == 1:
103        return 1, 1
104    elif n == 2:
105        return 2, 2
106    elif n == 3:
107        return 1, 2
108    elif n == 4:
109        return 0, 2
110    elif n == 5:
111        return 0, 1
112    else:
113        return None
114
115
116
117def index_to_voigt(i: int, j: int) -> Optional[int]:
118    """
119    Converts tensorial (2 indices) to Voigt's (1 index) notation
120
121    :param i: index in tensorial notation (0..2)
122    :type i: int
123    :param j: index in tensorial notation (0..2)
124    :type j: int
125    :return: index in Voigt's notation (0..6)
126    :rtype: int
127    """
128    return {v: k for k, v in VOIGT_INDEX_MAP.items()}.get((i, j))
129
130
131
132def transform_tensor(tensor: NDArray, arr: NDArray, fact_two: bool = False) -> NDArray:
133    """
134    Transforms tensor according to a given transformation matrix
135    :param tensor: tensor to be transformed
136    :type tensor: NDArray
137    :param arr: transformation matrix
138    :type arr: NDArray
139    :param fact_two: whether to apply factors 2, 4, ... during conversion to/from Voigt's notation
140    :type fact_two: bool
141    :return: transformed matrix
142    :rtype: NDArray
143    """
144    if tensor.shape in {(6, 6), (3, 3, 3, 3)}:
145        if tensor.shape == (6, 6):
146            voigt = True
147            tensor = from_voigt(tensor, fact_two)
148        else:
149            voigt = False
150        e = np.zeros(81).reshape((3, 3, 3 ,3))
151        for i, j, k, l in product(np.arange(3), np.arange(3), np.arange(3), np.arange(3)):
152            for a, b, c, d in product(np.arange(3), np.arange(3), np.arange(3), np.arange(3)):
153                e[i, j, k, l] += np.around(arr[i, a]*arr[j, b]*arr[k, c]*arr[l, d]*tensor[a, b, c, d], 2)
154        if voigt:
155            e = to_voigt(e, fact_two)
156        return e
157    else:
158        raise ValueError('Unknown shape of the input tensor.')
159
160
161
162def from_voigt(m: ndarray, div_two: bool = True) -> ndarray:
163    """
164    Transform the stress or strain array from Voigt's notation into a (3,3) matrix.
165    
166    In case it is a (6,) array, the function will produce a second-rank tensor (3x3 matrix).
167    In case it is a (6,6) matrix, the function will produce a fourth-rank tensor (3x3x3x3 matrix).
168
169    Parameters
170    ----------
171    m : ndarray
172        The stress or strain array in Voigt's notation.
173    div_two : bool, optional
174        Divide off-diagonal elements by two (default is True).
175
176    Returns
177    -------
178    ndarray
179        A transformed tensor.
180
181    Raises
182    ------
183    ValueError
184        If the input data has an unknown shape.
185    """
186
187    if m.shape == (6,):
188        # 2nd rank tensor: vector 6x1 -> matrix 3x3
189        e = m.copy()
190        if div_two:
191            e[3] /= 2.0
192            e[4] /= 2.0
193            e[5] /= 2.0
194        m = np.array([
195            [e[0], e[5], e[4]],
196            [e[5], e[1], e[3]],
197            [e[4], e[3], e[2]]
198        ])
199        return m
200    elif m.shape == (6, 6):
201        # 4th rank tensor: matrix 6x6 -> matrix 3x3x3x3
202        fact = np.ones(6)
203        if div_two:
204            fact[3:] = 0.5
205        e = np.zeros(81).reshape((3, 3, 3, 3))
206        for i, j in np.ndindex((6, 6)):
207            a, b = index_from_voigt(i)
208            c, d = index_from_voigt(j)
209            e[a, b, c, d] = m[i, j] * fact[i] * fact[j]
210            e[a, b, d, c] = m[i, j] * fact[i] * fact[j]
211            e[b, a, c, d] = m[i, j] * fact[i] * fact[j]
212            e[b, a, d, c] = m[i, j] * fact[i] * fact[j]
213        return e
214    else:
215        raise ValueError('Unknown shape of the input data')
216
217
218
219def to_voigt(m: NDArray, times_two: bool = True) -> NDArray:
220    """
221    Transform a stress or strain array tensor into Voigt's notation.
222    
223    For a second-rank tensor (3x3 matrix), it transforms into a Voigt's notation (6,).
224    For a fourth-rank tensor (3x3x3x3 matrix), it transforms into a Voigt's notation (6x6).
225
226    Parameters
227    ----------
228    m : NDArray
229        The stress or strain array tensor.
230    times_two : bool, optional
231        Multiply off-diagonal elements by two (default is True).
232
233    Returns
234    -------
235    NDArray
236        A transformed tensor.
237
238    Raises
239    ------
240    ValueError
241        If the input data has an unknown shape.
242    """
243    if m.shape == (3, 3):
244        voigt = np.array([
245            m[0, 0],
246            m[1, 1],
247            m[2, 2],
248            (m[2, 1] + m[1, 2]) / 2.0,
249            (m[2, 0] + m[0, 2]) / 2.0,
250            (m[1, 0] + m[0, 1]) / 2.0
251        ])
252        if times_two:
253            voigt[3] *= 2.0
254            voigt[4] *= 2.0
255            voigt[5] *= 2.0
256        return voigt
257    elif m.shape == (3, 3, 3, 3):
258        fact = np.ones(6)
259        if times_two:
260            fact[3:] = 2.0
261        voigt = np.zeros(36).reshape((6, 6))
262        for i, j in product(range(6), range(6)):
263            # index_from_voigt returns an Optional type, therefore if a "None" would occur 
264            # the function breaks at the unpacking step
265            a, b = index_from_voigt(i)
266            c, d = index_from_voigt(j)
267            voigt[i, j] = m[a, b, c, d] * fact[i] * fact[j]
268        return voigt
269    else:
270        raise ValueError('A stress matrix must be of shape (3,3)')
271
272
273def project_cubic(cij: NDArray) -> NDArray:
274    """
275    Computes to Cij tensor projected to cubic symmetry
276
277    :param cij: array of shape (6,6) the raw elasticity tensor
278    :type cij: NDArray
279    :return: the projected Cij tensor
280    :rtype: NDArray
281    """
282    cij = np.array(cij)
283    projected_cij = np.zeros((6, 6))
284    projected_cij[0, 0] = np.around(np.mean([cij[i, i] for i in np.arange(3)]), 2)
285    for i in np.arange(1, 3):
286        projected_cij[i, i] = projected_cij[0, 0]
287    projected_cij[3, 3] = np.around(np.mean([cij[i, i] for i in np.arange(3, 6)]), 2)
288    for i in np.arange(4, 6):
289        projected_cij[i, i] = projected_cij[3, 3]
290    projected_cij[1, 0] = np.around(np.mean([cij[1, 0], cij[2, 0], cij[2, 1]]), 2)
291    projected_cij[2, 0] = projected_cij[1, 0]
292    projected_cij[2, 1] = projected_cij[1, 0]
293    projected_cij[0, 1] = projected_cij[1, 0]
294    projected_cij[0, 2] = projected_cij[1, 0]
295    projected_cij[1, 2] = projected_cij[1, 0]
296    return projected_cij
297
298
299def project_hexagonal(cij: NDArray) -> NDArray:
300    """
301    Computes Cij tensor projected to hexagonal symmetry
302    :param cij: of shape (6,6) the raw elasticity tensor
303    :type cij: NDArray
304    :return: the projected Cij tensor
305    :rtype: NDArray
306    """
307
308    # convert from Cij to cij_hat, Moakher & Norris, Eq. 7
309    cij_hat = np.array(cij).copy()
310    for j in np.arange(3, 6):
311        for i in np.arange(3):
312            cij_hat[i, j] *= np.sqrt(2)
313            cij_hat[j, i] *= np.sqrt(2)
314        for i in np.arange(3, 6):
315            cij_hat[i, j] *= 2
316    # projecting on hexagonal symmetry, Moakher & Norris, Eq. A11a, A11b
317    c11st = (3*cij_hat[0,0]+3*cij_hat[1,1]+2*cij_hat[0,1]+2*cij_hat[5,5])/8.
318    c66st = (cij_hat[0,0]+cij_hat[1,1]-2*cij_hat[0,1]+2*cij_hat[5,5])/4.
319    # Moakher & Norris, Eq. A15         
320    projected_cij = np.zeros((6, 6))
321    projected_cij[0, 0] = np.around(c11st, 2)
322    projected_cij[1, 1] = projected_cij[0, 0]
323    projected_cij[2, 2] = np.around(cij_hat[2, 2], 2)
324    projected_cij[3, 3] = np.around(np.mean([cij_hat[3, 3], cij_hat[4, 4]]), 2)
325    projected_cij[4, 4] = projected_cij[3, 3]
326    projected_cij[5, 5] = np.around(c66st, 2)
327    projected_cij[1, 0] = np.around(c11st-c66st, 2)
328    projected_cij[0, 1] = projected_cij[1, 0]
329    projected_cij[2, 0] = np.around(np.mean([cij_hat[2, 0], cij_hat[2, 1]]), 2)
330    projected_cij[2, 1] = projected_cij[2, 0]
331    projected_cij[0, 2] = projected_cij[2, 0]
332    projected_cij[1, 2] = projected_cij[2, 0]
333    # convert back,Moakher & Norris, Eq. 7
334    for j in range(3, 6):
335        for i in range(3):
336            projected_cij[i, j] /= np.sqrt(2)
337            projected_cij[j, i] /= np.sqrt(2)
338        for i in range(3, 6):
339            projected_cij[i, j] /= 2
340
341    return projected_cij
342
343
344
345###############################################################################################################
346# LEFT OVER CODE FROM THE PAST, PERHAPS REIMPLEMENT IN SOME FORM??? ---> move to scopes.cij?
347    
348#     def read_data_from_VASP_Cij_calculation(self, case,
349#                                             path_to_CONTCAR = ['V=eq', 'CONTCAR'],
350#                                             path_to_Cij = ['Cij', 'eps=0.014', 'Cij_cubic.dat']):
351#         """
352#         Reads and processes the data needed for Debye model.
353
354#         Arguments:
355#         case   path to the folder containing <eq> and <Cij/eps=0.014> folders
356#         """
357
358#         from pymatgen.io.vasp import Poscar, Oszicar
359#         from pymatgen.core.periodic_table import Element
360#         from os.path import join
361#         from numpy import loadtxt
362#         from numpy.linalg import inv
363
364#         struct = Poscar.from_file(join(case, *path_to_CONTCAR)).structure
365#         self.data['V'] = struct.volume/struct.num_sites
366#         E0 = Oszicar(join(case, *path_to_CONTCAR[:-1], 'OSZICAR')).final_energy
367#         self.data['E0'] = E0/struct.num_sites
368#         self.struct = struct
369        
370#         Cij = loadtxt(join(case, *path_to_Cij))
371#         Sij = inv(Cij)
372# #         https://wiki.materialsproject.org/Elasticity_calculations
373#         BV = sum([sum([Cij[i][j] for j in range(3)]) for i in range(3)])/9
374#         BR = 1/sum([sum([Sij[i][j] for j in range(3)]) for i in range(3)])
375#         BH = 0.5*(BV+BR)
376#         GV = (sum(Cij[i][i] for i in range(3)) + sum(3*Cij[i][i] for i in range(3, 6)) - (Cij[0][1]+Cij[1][2]+Cij[0][2]))/15
377#         GR = 15/(4*sum(Sij[i][i] for i in range(3)) + sum(3*Sij[i][i] for i in range(3, 6)) - 4*(Sij[0][1]+Sij[1][2]+Sij[0][2]))
378#         GH = 0.5*(GV+GR)
379#         self.data['B'] = BH
380#         self.data['nu'] = (3*BH-2*GH)/(6*BH+2*GH)
381        
382            
383#         M = sum([Element[s.species_string].data['Atomic mass'] for s in struct.sites])
384#         # atomic mass unit -> g:  u = 1.66054e-24 g
385#         # volume: Ang^3 = 1e-30 m^3 = 1e-24 cm^3
386#         self.data['rho'] = 1.66054*M/struct.volume
387###############################################################################################################