1# FHIaims.py - IO routines for phonopy-FHI-aims
2# methods compatible with the corresponding ones from ase.io.aims
3# only minimal subset of functionality required within phonopy context is implemented
4#
5# Copyright (C) 2009-2011 Joerg Meyer (jm)
6# All rights reserved.
7#
8# This file is part of phonopy.
9#
10# Redistribution and use in source and binary forms, with or without
11# modification, are permitted provided that the following conditions
12# are met:
13#
14# * Redistributions of source code must retain the above copyright
15#   notice, this list of conditions and the following disclaimer.
16#
17# * Redistributions in binary form must reproduce the above copyright
18#   notice, this list of conditions and the following disclaimer in
19#   the documentation and/or other materials provided with the
20#   distribution.
21#
22# * Neither the name of the phonopy project nor the names of its
23#   contributors may be used to endorse or promote products derived
24#   from this software without specific prior written permission.
25#
26# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
29# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
30# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
31# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
32# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
33# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
34# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
35# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
36# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
37# POSSIBILITY OF SUCH DAMAGE.
38#
39# Modified 2020 by Florian Knoop
40
41import sys
42import numpy as np
43from phonopy.structure.atoms import PhonopyAtoms as Atoms
44from phonopy.interface.vasp import check_forces, get_drift_forces
45
46
47# FK 2018/07/19
48def lmap(func, lis):
49    """Python2/3 compatibility:
50        replace map(int, list) with lmap(int, list) that always returns a list
51        instead of an iterator. Otherwise conflicts with np.array in python3. """
52    return list(map(func, lis))
53
54
55def read_aims(filename):
56    """Method to read FHI-aims geometry files in phonopy context."""
57
58    lines = open(filename, "r").readlines()
59
60    cell = []
61    is_frac = []
62    positions = []
63    symbols = []
64    magmoms = []
65    for line in lines:
66        fields = line.split()
67        if not len(fields):
68            continue
69        if fields[0] == "lattice_vector":
70            vec = lmap(float, fields[1:4])
71            cell.append(vec)
72        elif fields[0][0:4] == "atom":
73            if fields[0] == "atom":
74                frac = False
75            elif fields[0] == "atom_frac":
76                frac = True
77            pos = lmap(float, fields[1:4])
78            sym = fields[4]
79            is_frac.append(frac)
80            positions.append(pos)
81            symbols.append(sym)
82            magmoms.append(None)
83        # implicitly assuming that initial_moments line adhere to FHI-aims geometry.in specification,
84        # i.e. two subsequent initial_moments lines do not occur
85        # if they do, the value specified in the last line is taken here - without any warning
86        elif fields[0] == "initial_moment":
87            magmoms[-1] = float(fields[1])
88
89    for (n, frac) in enumerate(is_frac):
90        if frac:
91            pos = [
92                sum([positions[n][l] * cell[l][i] for l in range(3)]) for i in range(3)
93            ]
94            positions[n] = pos
95    if None in magmoms:
96        atoms = Atoms(cell=cell, symbols=symbols, positions=positions)
97    else:
98        atoms = Atoms(cell=cell, symbols=symbols, positions=positions, magmoms=magmoms)
99
100    return atoms
101
102
103def write_aims(filename, atoms):
104    """Method to write FHI-aims geometry files in phonopy context."""
105
106    lines = ""
107    lines += "# geometry.in for FHI-aims \n"
108    lines += "# | generated by phonopy.FHIaims.write_aims() \n"
109
110    lattice_vector_line = "lattice_vector " + "%16.16f " * 3 + "\n"
111    for vec in atoms.get_cell():
112        lines += lattice_vector_line % tuple(vec)
113
114    N = atoms.get_number_of_atoms()
115
116    atom_line = "atom " + "%16.16f " * 3 + "%s \n"
117    positions = atoms.get_positions()
118    symbols = atoms.get_chemical_symbols()
119
120    initial_moment_line = "initial_moment %16.6f\n"
121    magmoms = atoms.get_magnetic_moments()
122
123    for n in range(N):
124        lines += atom_line % (tuple(positions[n]) + (symbols[n],))
125        if magmoms is not None:
126            lines += initial_moment_line % magmoms[n]
127
128    with open(filename, "w") as f:
129        f.write(lines)
130
131
132class Atoms_with_forces(Atoms):
133    """ Hack to phonopy.atoms to maintain ASE compatibility also for forces."""
134
135    def get_forces(self):
136        return self.forces
137
138
139def read_aims_output(filename):
140    """ Read FHI-aims output and
141        return geometry, energy and forces from last self-consistency iteration"""
142
143    lines = open(filename, "r").readlines()
144
145    l = 0
146    N = 0
147    while l < len(lines):
148        line = lines[l]
149        if "| Number of atoms" in line:
150            N = int(line.split()[5])
151        elif "| Unit cell:" in line:
152            cell = []
153            for i in range(3):
154                l += 1
155                vec = lmap(float, lines[l].split()[1:4])
156                cell.append(vec)
157        elif ("Atomic structure:" in line) or ("Updated atomic structure:" in line):
158            if "Atomic structure:" in line:
159                i_sym = 3
160                i_pos_min = 4
161                i_pos_max = 7
162            elif "Updated atomic structure:" in line:
163                i_sym = 4
164                i_pos_min = 1
165                i_pos_max = 4
166            l += 1
167            symbols = []
168            positions = []
169            for n in range(N):
170                l += 1
171                fields = lines[l].split()
172                sym = fields[i_sym]
173                pos = lmap(float, fields[i_pos_min:i_pos_max])
174                symbols.append(sym)
175                positions.append(pos)
176        elif "Total atomic forces" in line:
177            forces = []
178            for i in range(N):
179                l += 1
180                force = lmap(float, lines[l].split()[-3:])
181                forces.append(force)
182        l += 1
183
184    atoms = Atoms_with_forces(cell=cell, symbols=symbols, positions=positions)
185    atoms.forces = forces
186
187    return atoms
188
189
190def write_supercells_with_displacements(supercell,
191                                        cells_with_disps,
192                                        ids,
193                                        pre_filename="geometry.in",
194                                        width=3):
195    """Writes perfect supercell and supercells with displacements
196
197    Args:
198        supercell: perfect supercell
199        cells_with_disps: supercells with displaced atoms
200        filename: root-filename
201    """
202
203    # original cell
204    write_aims(pre_filename + ".supercell", supercell)
205
206    # displaced cells
207    for i, cell in zip(ids, cells_with_disps):
208        filename = "{pre_filename}-{0:0{width}}".format(
209            i, pre_filename=pre_filename, width=width)
210        write_aims(filename, cell)
211
212
213def parse_set_of_forces(num_atoms, forces_filenames, verbose=True):
214    """parse the forces from output files in `forces_filenames`"""
215    is_parsed = True
216    force_sets = []
217    for i, filename in enumerate(forces_filenames):
218        if verbose:
219            sys.stdout.write("%d. " % (i + 1))
220
221        atoms = read_aims_output(filename)
222        forces = atoms.forces
223        if check_forces(forces, num_atoms, filename, verbose=verbose):
224            drift_force = get_drift_forces(forces, filename=filename, verbose=verbose)
225            force_sets.append(np.array(forces) - drift_force)
226        else:
227            is_parsed = False
228
229    if is_parsed:
230        return force_sets
231    else:
232        return []
233