1# Copyright (C) 2011 Atsushi Togo
2# All rights reserved.
3#
4# This file is part of phonopy.
5#
6# Redistribution and use in source and binary forms, with or without
7# modification, are permitted provided that the following conditions
8# are met:
9#
10# * Redistributions of source code must retain the above copyright
11#   notice, this list of conditions and the following disclaimer.
12#
13# * Redistributions in binary form must reproduce the above copyright
14#   notice, this list of conditions and the following disclaimer in
15#   the documentation and/or other materials provided with the
16#   distribution.
17#
18# * Neither the name of the phonopy project nor the names of its
19#   contributors may be used to endorse or promote products derived
20#   from this software without specific prior written permission.
21#
22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
33# POSSIBILITY OF SUCH DAMAGE.
34
35import numpy as np
36
37directions_axis = np.array([[1, 0, 0],
38                            [0, 1, 0],
39                            [0, 0, 1]])
40
41directions_diag = np.array([[1, 0, 0],
42                            [0, 1, 0],
43                            [0, 0, 1],
44                            [1, 1, 0],
45                            [1, 0, 1],
46                            [0, 1, 1],
47                            [1, -1, 0],
48                            [1, 0, -1],
49                            [0, 1, -1],
50                            [1, 1, 1],
51                            [1, 1, -1],
52                            [1, -1, 1],
53                            [-1, 1, 1]])
54
55
56def directions_to_displacement_dataset(displacement_directions,
57                                       distance,
58                                       supercell):
59    lattice = supercell.get_cell()
60    first_atoms = []
61    for disp in displacement_directions:
62        direction = disp[1:]
63        disp_cartesian = np.dot(direction, lattice)
64        disp_cartesian *= distance / np.linalg.norm(disp_cartesian)
65        first_atoms.append({'number': int(disp[0]),
66                            'displacement': disp_cartesian.tolist()})
67    displacement_dataset = {
68        'natom': supercell.get_number_of_atoms(),
69        'first_atoms': first_atoms}
70
71    return displacement_dataset
72
73
74def get_least_displacements(symmetry,
75                            is_plusminus='auto',
76                            is_diagonal=True,
77                            is_trigonal=False,
78                            log_level=0):
79    """Return a set of displacements
80
81    Returns
82    -------
83    array_like
84        List of directions with respect to axes. This gives only the
85        symmetrically non equivalent directions. The format is like:
86           [[0, 1, 0, 0],
87            [7, 1, 0, 1], ...]
88        where each list is defined by:
89           First value:      Atom index in supercell starting with 0
90           Second to fourth: If the direction is displaced or not (1, 0, or -1)
91                             with respect to the axes.
92
93    """
94    displacements = []
95    if is_diagonal:
96        directions = directions_diag
97    else:
98        directions = directions_axis
99
100    if log_level > 2:
101        print("Site point symmetry:")
102
103    for atom_num in symmetry.get_independent_atoms():
104        site_symmetry = symmetry.get_site_symmetry(atom_num)
105
106        if log_level > 2:
107            print("Atom %d" % (atom_num + 1))
108            for i, rot in enumerate(site_symmetry):
109                print("----%d----" % (i + 1))
110                for v in rot:
111                    print("%2d %2d %2d" % tuple(v))
112
113        for disp in get_displacement(site_symmetry,
114                                     directions,
115                                     is_trigonal,
116                                     log_level):
117            displacements.append([atom_num,
118                                  disp[0], disp[1], disp[2]])
119            if is_plusminus == 'auto':
120                if is_minus_displacement(disp, site_symmetry):
121                    displacements.append([atom_num,
122                                          -disp[0], -disp[1], -disp[2]])
123            elif is_plusminus is True:
124                displacements.append([atom_num,
125                                      -disp[0], -disp[1], -disp[2]])
126
127    return displacements
128
129
130def get_displacement(site_symmetry,
131                     directions=directions_diag,
132                     is_trigonal=False,
133                     log_level=0):
134    # One
135    sitesym_num, disp = get_displacement_one(site_symmetry,
136                                             directions)
137    if disp is not None:
138        if log_level > 2:
139            print("Site symmetry used to expand a direction %s" % disp[0])
140            print(site_symmetry[sitesym_num])
141        return disp
142    # Two
143    sitesym_num, disps = get_displacement_two(site_symmetry,
144                                              directions)
145    if disps is not None:
146        if log_level > 2:
147            print("Site symmetry used to expand directions %s %s" %
148                  (disps[0], disps[1]))
149            print(site_symmetry[sitesym_num])
150
151        if is_trigonal:
152            disps_new = [disps[0]]
153            if is_trigonal_axis(site_symmetry[sitesym_num]):
154                if log_level > 2:
155                    print("Trigonal axis is found.")
156                disps_new.append(np.dot(disps[0],
157                                        site_symmetry[sitesym_num].T))
158                disps_new.append(np.dot(disps_new[1],
159                                        site_symmetry[sitesym_num].T))
160            disps_new.append(disps[1])
161            return disps_new
162        else:
163            return disps
164    # Three
165    return [directions[0], directions[1], directions[2]]
166
167
168def get_displacement_one(site_symmetry,
169                         directions=directions_diag):
170    for direction in directions:
171        rot_directions = []
172        for r in site_symmetry:
173            rot_directions.append(np.dot(direction, r.T))
174        num_sitesym = len(site_symmetry)
175        for i in range(num_sitesym):
176            for j in range(i+1, num_sitesym):
177                det = determinant(direction,
178                                  rot_directions[i],
179                                  rot_directions[j])
180                if det != 0:
181                    return i, [direction]
182    return None, None
183
184
185def get_displacement_two(site_symmetry,
186                         directions=directions_diag):
187    for direction in directions:
188        rot_directions = []
189        for r in site_symmetry:
190            rot_directions.append(np.dot(direction, r.T))
191        num_sitesym = len(site_symmetry)
192        for i in range(num_sitesym):
193            for second_direction in directions:
194                det = determinant(direction,
195                                  rot_directions[i],
196                                  second_direction)
197                if det != 0:
198                    return i, [direction, second_direction]
199    return None, None
200
201
202def is_minus_displacement(direction, site_symmetry):
203    is_minus = True
204    for r in site_symmetry:
205        rot_direction = np.dot(direction, r.T)
206        if (rot_direction + direction).any():
207            continue
208        else:
209            is_minus = False
210            break
211    return is_minus
212
213
214def is_trigonal_axis(r):
215    r3 = np.dot(np.dot(r, r), r)
216    if (r3 == np.eye(3, dtype=int)).all():
217        return True
218    else:
219        return False
220
221
222def determinant(a, b, c):
223    det = (a[0] * b[1] * c[2] - a[0] * b[2] * c[1]
224           + a[1] * b[2] * c[0] - a[1] * b[0] * c[2]
225           + a[2] * b[0] * c[1] - a[2] * b[1] * c[0])
226    return det
227
228
229def get_random_displacements_dataset(num_supercells,
230                                     distance,
231                                     num_atoms,
232                                     random_seed=None):
233    if (np.issubdtype(type(random_seed), np.integer) and
234        random_seed >= 0 and random_seed < 2 ** 32):
235        seed = random_seed
236    else:
237        seed = None
238
239    disps = get_random_directions(num_atoms * num_supercells,
240                                  random_seed=random_seed) * distance
241    supercell_disps = np.array(disps.reshape(num_supercells, num_atoms, 3),
242                               dtype='double', order='C')
243    dataset = {'displacements': supercell_disps}
244
245    if seed is not None:
246        dataset['random_seed'] = seed
247    return dataset
248
249
250def get_random_directions(num_atoms, random_seed=None):
251    """Returns random directions in sphere with radius 1"""
252
253    if (np.issubdtype(type(random_seed), np.integer) and
254        random_seed >= 0 and random_seed < 2 ** 32):
255        np.random.seed(random_seed)
256
257    xy = np.random.randn(3, num_atoms)
258    r = np.sqrt((xy ** 2).sum(axis=0))
259    return (xy / r).T
260
261
262def print_displacements(symmetry,
263                        directions=directions_diag):
264    displacements = get_least_displacements(symmetry, directions)
265    print("Least displacements:")
266    print("Atom       Directions")
267    print("----------------------------")
268    for key in displacements:
269        print("%4d  %s" % (key + 1, displacements[key]))
270