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