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 sys
36import numpy as np
37from phonopy.interface.vasp import get_drift_forces, check_forces
38from phonopy.structure.atoms import PhonopyAtoms as Atoms
39from phonopy.structure.symmetry import Symmetry
40from phonopy.structure.cells import get_angles, get_cell_parameters
41from phonopy.harmonic.force_constants import similarity_transformation
42
43
44def parse_set_of_forces(disps,
45                        forces_filenames,
46                        supercell,
47                        wien2k_P1_mode=False,  # Only for the test
48                        symmetry_tolerance=None,
49                        verbose=True):
50    if symmetry_tolerance is None:
51        symprec = 1e-5
52    else:
53        symprec = symmetry_tolerance
54
55    num_atoms = supercell.get_number_of_atoms()
56    lattice = supercell.get_cell()
57    is_parsed = True
58    force_sets = []
59
60    for i, (filename, disp) in enumerate(zip(forces_filenames, disps)):
61        if verbose:
62            sys.stdout.write("%d. " % (i + 1))
63
64        # Parse wien2k case.scf file
65        wien2k_forces = _get_forces_wien2k(filename, lattice)
66        if not wien2k_P1_mode:
67            forces = _distribute_forces(
68                supercell,
69                disp,
70                wien2k_forces,
71                filename,
72                symprec)
73        else:
74            if num_atoms != len(wien2k_forces):
75                forces = []
76            else:
77                forces = wien2k_forces
78
79        if check_forces(forces, num_atoms, filename, verbose=verbose):
80            drift_force = get_drift_forces(forces,
81                                           filename=filename,
82                                           verbose=verbose)
83            force_sets.append(np.array(forces) - drift_force)
84        else:
85            is_parsed = False
86
87    if is_parsed:
88        return force_sets
89    else:
90        return []
91
92
93def parse_wien2k_struct(filename):
94    with open(filename) as f:
95        # 1
96        title = f.readline().rstrip()
97
98        # 2
99        num_site = int(f.readline()[27:30])
100
101        # 3
102        f.readline()
103
104        # 4
105        line = f.readline()
106        a = float(line[0:10])
107        b = float(line[10:20])
108        c = float(line[20:30])
109        alpha = float(line[30:40])
110        beta = float(line[40:50])
111        gamma = float(line[50:60])
112
113        lattice = _transform_axis(alpha, beta, gamma, a, b, c)
114
115        symbols = []
116        positions = []
117        npts = []
118        r0s = []
119        rmts = []
120
121        for i in range(num_site):
122            # 5
123            line = f.readline()
124            x = float(line[12:22])
125            y = float(line[25:35])
126            z = float(line[38:48])
127            positions.append([x, y, z])
128
129            # 6
130            line = f.readline()
131            multi = int(line[15:17])
132
133            for j in range(multi - 1):
134                line = f.readline()
135                x = float(line[12:22])
136                y = float(line[25:35])
137                z = float(line[38:48])
138                positions.append([x, y, z])
139
140            # 7
141            line = f.readline()
142            chemical_symbol = line[0:2].strip()
143            npt = int(line[15:20])
144            r0 = float(line[25:35])
145            rmt = float(line[40:50])
146
147            for j in range(multi):
148                symbols.append(chemical_symbol)
149                npts.append(npt)
150                r0s.append(r0)
151                rmts.append(rmt)
152
153            # 8 - 10
154            for j in range(3):
155                f.readline()
156
157        cell = Atoms(symbols=symbols,
158                     scaled_positions=positions,
159                     cell=lattice)
160
161        return cell, npts, r0s, rmts
162
163
164def write_supercells_with_displacements(supercell,
165                                        cells_with_displacements,
166                                        ids,
167                                        npts, r0s, rmts,
168                                        num_unitcells_in_supercell,
169                                        pre_filename="wien2k",
170                                        width=3):
171    npts_super = []
172    r0s_super = []
173    rmts_super = []
174    for i, j, k in zip(npts, r0s, rmts):
175        for l in range(num_unitcells_in_supercell):
176            npts_super.append(i)
177            r0s_super.append(j)
178            rmts_super.append(k)
179
180    _pre_filename = pre_filename.split('/')[-1] + "S"
181    write_wein2k(_pre_filename,
182                 supercell,
183                 npts_super,
184                 r0s_super,
185                 rmts_super)
186
187    for i, cell in zip(ids, cells_with_displacements):
188        symmetry = Symmetry(cell)
189        filename = "{pre_filename}-{0:0{width}}.in".format(
190            i, pre_filename=_pre_filename, width=width)
191        print("Number of non-equivalent atoms in %s: %d" %
192              (filename, len(symmetry.get_independent_atoms())))
193        write_wein2k(filename,
194                     cell,
195                     npts_super,
196                     r0s_super,
197                     rmts_super)
198
199
200def write_wein2k(filename, cell, npts, r0s, rmts):
201    with open(filename, 'w') as w:
202        w.write(_get_wien2k_struct(cell, npts, r0s, rmts))
203
204
205def _get_wien2k_struct(cell, npts, r0s, rmts):
206    num_atom = cell.get_number_of_atoms()
207    lattice = cell.get_cell()
208    a, b, c = get_cell_parameters(lattice)
209    alpha, beta, gamma = get_angles(lattice)
210    positions = cell.get_scaled_positions()
211    symbols = cell.get_chemical_symbols()
212    numbers = cell.get_atomic_numbers()
213
214    text = ""
215
216    # 1
217    text += "Title\n"
218
219    # 2
220    text += "%-4s%23s%3d\n" % ("P", "LATTICE,NONEQUIV.ATOMS:", num_atom)
221
222    # 3
223    text += "%13s%4s\n" % ("MODE OF CALC=", "RELA")
224
225    # 4
226    text += "%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f\n" % (
227        a, b, c, alpha, beta, gamma)
228
229    for i, pos in enumerate(positions):
230        for j in (0, 1, 2):
231            if pos[j] < 0:
232                pos[j] += 1
233            if int(float("%10.8f" % pos[j])) == 1:
234                pos[j] = 0.0
235
236        # 5 format (4X,I4,4X,F10.8,3X,F10.8,3X,F10.8)
237        text += "%4s%4d%4s%10.8f%3s%10.8f%3s%10.8f\n" % (
238            "ATOM", -(i + 1), ": X=", pos[0], " Y=", pos[1], " Z=", pos[2])
239
240        # 6  format (15X,I2,17X,I2)
241        text += "%15s%2d%17s%2d\n" % ("MULT=", 1, "ISPLIT=", 8)
242
243        # 7 format (A10,5X,I5,5X,F10.8,5X,F10.5,5X,F5.2)
244        npt = npts[i]
245        r0 = r0s[i]
246        rmt = rmts[i]
247        text += "%-10s%5s%5d%5s%10.8f%5s%10.5f%5s%5.1f\n" % (
248            symbols[i], "NPT=", npt, "R0=", r0, "RMT=", rmt, "Z:", numbers[i])
249
250        # 8 - 10 format (20X,3F10.7)
251        text += "%-20s%10.7f%10.7f%10.7f\n" % ("LOCAL ROT MATRIX:", 1, 0, 0)
252        text += "%-20s%10.7f%10.7f%10.7f\n" % ("", 0, 1, 0)
253        text += "%-20s%10.7f%10.7f%10.7f\n" % ("", 0, 0, 1)
254
255    text += "   0      NUMBER OF SYMMETRY OPERATIONS"
256
257    return text
258
259
260def _transform_axis(alpha, beta, gamma, a, b, c):
261    alpha = alpha / 180 * np.pi
262    beta = beta / 180 * np.pi
263    gamma = gamma / 180 * np.pi
264
265    cz = c
266
267    bz = np.cos(alpha) * b
268    by = np.sin(alpha) * b
269
270    az = np.cos(beta) * a
271    ay = (np.cos(gamma) - np.cos(beta) * np.cos(alpha)) / np.sin(alpha) * a
272    ax = np.sqrt(a**2 - ay**2 - az**2)
273
274    return [ax, ay, az], [0, by, bz], [0, 0, cz]
275
276
277def _parse_core_param(file):
278    npts = []
279    r0s = []
280    rmts = []
281
282    for line in file:
283        if line.strip()[0] == '#':
284            continue
285
286        vals = line.strip().split()
287        npts.append(int(vals[1]))
288        r0s.append(float(vals[2]))
289        rmts.append(float(vals[3]))
290
291    return npts, r0s, rmts
292
293
294def _get_forces_wien2k(filename, lattice):
295    forces = []
296    red_lattice = []
297
298    for v in lattice:
299        red_lattice.append(v / np.sqrt(np.vdot(v, v)))
300
301    num_atom = 0
302    for line in open(filename):
303        if line.count('total forces') > 0:
304            if line[:4] == ":FGL":
305                fx = float(line[29:45])
306                fy = float(line[45:61])
307                fz = float(line[61:77])
308                forces.append(np.dot([fx, fy, fz], red_lattice))
309                num_atom = int(line[4:7])
310
311    return forces[-num_atom:]
312
313
314def _distribute_forces(supercell, disp, forces, filename, symprec):
315    natom = supercell.get_number_of_atoms()
316    lattice = supercell.get_cell()
317    symbols = supercell.get_chemical_symbols()
318    positions = supercell.get_positions() + disp
319    cell = Atoms(cell=lattice,
320                 positions=positions,
321                 symbols=symbols,
322                 pbc=True)
323    symmetry = Symmetry(cell, symprec)
324    independent_atoms = symmetry.get_independent_atoms()
325
326    # Rotation matrices in Cartesian
327    rotations = []
328    for r in symmetry.get_symmetry_operations()['rotations']:
329        rotations.append(similarity_transformation(lattice.T, r))
330
331    map_operations = symmetry.get_map_operations()
332    map_atoms = symmetry.get_map_atoms()
333
334    atoms_in_dot_scf = _get_independent_atoms_in_dot_scf(filename)
335
336    if len(forces) != len(atoms_in_dot_scf):
337        print("%s does not contain necessary information." % filename)
338        print("Plese check if there are \"FGL\" lines with")
339        print("\"total forces\" are required.")
340        return False
341
342    if len(atoms_in_dot_scf) == natom:
343        print("It is assumed that there is no symmetrically-equivalent "
344              "atoms in ")
345        print("\'%s\' at wien2k calculation." % filename)
346        force_set = forces
347    elif len(forces) != len(independent_atoms):
348        print("Non-equivalent atoms of %s could not be recognized by phonopy."
349              % filename)
350        return False
351    else:
352        # 1. Transform wien2k forces to those on independent atoms
353        indep_atoms_to_wien2k = []
354        forces_remap = []
355        for i, pos_wien2k in enumerate(atoms_in_dot_scf):
356            for j, pos in enumerate(cell.get_scaled_positions()):
357                diff = pos_wien2k - pos
358                diff -= np.rint(diff)
359                if (abs(diff) < symprec).all():
360                    forces_remap.append(
361                        np.dot(rotations[map_operations[j]], forces[i]))
362                    indep_atoms_to_wien2k.append(map_atoms[j])
363                    break
364
365        if len(forces_remap) != len(forces):
366            print("Atomic position mapping between Wien2k and phonopy failed.")
367            print("If you think this is caused by a bug of phonopy")
368            print("please report it in the phonopy mainling list.")
369            return False
370
371        # 2. Distribute forces from independent to dependent atoms.
372        force_set = []
373        for i in range(natom):
374            j = indep_atoms_to_wien2k.index(map_atoms[i])
375            force_set.append(np.dot(
376                rotations[map_operations[i]].T, forces_remap[j]))
377
378    return force_set
379
380
381def _get_independent_atoms_in_dot_scf(filename):
382    positions = []
383    for line in open(filename):
384        if line[:4] == ":POS":
385            if "POSITION" in line:
386                x = float(line[30:37])
387                y = float(line[38:45])
388                z = float(line[46:53])
389            else:
390                x = float(line[27:34])
391                y = float(line[35:42])
392                z = float(line[43:50])
393            num_atom = int(line[4:7])
394            positions.append([x, y, z])
395
396    return np.array(positions)[-num_atom:]
397
398
399if __name__ == '__main__':
400    from optparse import OptionParser
401    from phonopy.interface.vasp import write_vasp, read_vasp
402
403    def clean_scaled_positions(cell):
404        positions = cell.get_scaled_positions()
405        for pos in positions:
406            for i in (0, 1, 2):
407                # The following %19.16f follows write_vasp
408                if float("%19.16f" % pos[i]) >= 1:
409                    pos[i] -= 1.0
410        cell.set_scaled_positions(positions)
411
412    parser = OptionParser()
413    parser.set_defaults(w2v=False, v2w=False)
414    parser.add_option("-w", dest="w2v",
415                      action="store_true",
416                      help="Convert WIEN2k to VASP")
417    parser.add_option("-v", dest="v2w",
418                      action="store_true",
419                      help="Convert VASP to WIEN2k")
420    (options, args) = parser.parse_args()
421
422    from phonopy.units import Bohr
423
424    if options.v2w:
425        cell = read_vasp(args[0])
426        lattice = cell.get_cell() / Bohr
427        cell.set_cell(lattice)
428        npts, r0s, rmts = _parse_core_param(open(args[1]))
429        text = _get_wien2k_struct(cell, npts, r0s, rmts)
430        print(text)
431
432    elif options.w2v:
433        cell, npts, r0s, rmts = parse_wien2k_struct(args[0])
434        positions = cell.get_scaled_positions()
435        lattice = cell.get_cell() * Bohr
436        cell.set_cell(lattice)
437        cell.set_scaled_positions(positions)
438        clean_scaled_positions(cell)
439        write_vasp("POSCAR.wien2k", cell, direct=True)
440        w = open("wien2k_core.dat", 'w')
441
442        w.write("# symbol       npt       r0             rmt\n")
443        for symbol, npt, r0, rmt in \
444                zip(cell.get_chemical_symbols(), npts, r0s, rmts):
445            w.write("%-10s     %5d     %10.8f     %10.5f\n" %
446                    (symbol, npt, r0, rmt))
447    else:
448        print("You need to set -r or -w option.")
449