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