1# Copyright (C) 2019 Antti J. Karttunen (antti.j.karttunen@iki.fi) 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 os 37import numpy as np 38 39from phonopy.interface.vasp import check_forces, get_drift_forces 40from phonopy.structure.atoms import PhonopyAtoms as Atoms 41 42 43def parse_set_of_forces(num_atoms, 44 forces_filenames, 45 verbose=True): 46 # Filenames = subdirectories supercell-001, supercell-002, ... 47 force_sets = [] 48 for i, filename in enumerate(forces_filenames): 49 if verbose: 50 sys.stdout.write("%d. " % (i + 1)) 51 52 f_gradient = open(os.path.join(filename, 'gradient')) 53 lines = f_gradient.readlines() 54 f_gradient.close() 55 # Structure of the gradient file: 56 #$grad cartesian gradients 57 # cycle = 1 SCF energy = -578.5931883878 |dE/dxyz| = 0.000007 58 # coordinates (num_atoms lines) 59 # gradients (num_atoms lines) 60 #$end 61 turbomole_forces = [] 62 for line in lines[2 + num_atoms : 2 + 2 * num_atoms]: 63 # Replace D with E in double precision floats 64 turbomole_forces.append([float(x.replace('D', 'E')) for x in line.split()]) 65 66 # Change from gradient to force by inverting the sign 67 # Units: hartree / Bohr 68 turbomole_forces = np.negative(turbomole_forces) 69 70 if check_forces(turbomole_forces, num_atoms, filename, verbose=verbose): 71 is_parsed = True 72 drift_force = get_drift_forces(turbomole_forces, 73 filename=filename, 74 verbose=verbose) 75 force_sets.append(np.array(turbomole_forces) - drift_force) 76 else: 77 is_parsed = False 78 79 if is_parsed: 80 return force_sets 81 else: 82 return [] 83 84 85def read_turbomole(filename): 86 # filename is typically "control" 87 f_turbomole = open(filename) 88 turbomole_in = TurbomoleIn(f_turbomole.readlines()) 89 f_turbomole.close() 90 tags = turbomole_in.get_tags() 91 cell = Atoms(cell=tags['lattice_vectors'], 92 symbols=tags['atomic_species'], 93 positions=tags['coordinates']) 94 95 return cell 96 97 98def write_turbomole(filename, cell): 99 # Write geometry in a new directory 100 # Check if directory exists (directory supercell will already exist for phono3py) 101 if not os.path.exists(filename): 102 os.mkdir(filename) 103 104 # Create control file 105 lines = '$title ' + filename + '\n' 106 lines += '$symmetry c1\n' 107 lines += '$coord file=coord\n' 108 lines += '$periodic 3\n' 109 lines += '$kpoints\n' 110 lines += ' nkpoints KPOINTS_HERE\n' 111 lines += '$scfconv 10\n' 112 lines += '$lattice\n' 113 lattice = cell.get_cell() 114 for lattvec in lattice: 115 lines += ("%12.8f" * 3 + "\n") % tuple(lattvec) 116 lines += '$end\n' 117 f_control = open(os.path.join(filename, 'control'), 'w') 118 f_control.write(lines) 119 f_control.close() 120 121 # Create coord file 122 symbols = cell.get_chemical_symbols() 123 positions = cell.get_positions() 124 lines = '$coord\n' 125 for atom, pos in zip(symbols, positions): 126 lines += ("%16.12f"*3 + " %s\n") % (pos[0], pos[1], pos[2], atom.lower()) 127 lines += '$end\n' 128 f_coord = open(os.path.join(filename, 'coord'), 'w') 129 f_coord.write(lines) 130 f_coord.close() 131 132 133def write_supercells_with_displacements(supercell, 134 cells_with_displacements, 135 ids, 136 pre_filename="supercell", 137 width=3): 138 139 write_turbomole(pre_filename, supercell) 140 for i, cell in zip(ids, cells_with_displacements): 141 filename = "{pre_filename}-{0:0{width}}".format( 142 i, pre_filename=pre_filename, width=width) 143 write_turbomole(filename, cell) 144 145 146class TurbomoleIn: 147 def __init__(self, lines): 148 self._tags = {'lattice_vectors': None, 149 'atomic_species': None, 150 'coordinates': None} 151 152 self._values = None 153 self._collect(lines) 154 155 def get_tags(self): 156 return self._tags 157 158 def _collect(self, lines): 159 # Reads TURBOMOLE control and coord files (lattice vectors, cartesian atomic positions). 160 l = 0 161 lattvecs = [] 162 aspecies = [] 163 coords = [] 164 while l < len(lines): 165 line = lines[l] 166 # Look for lattice vectors 167 # Only supports atomic units, $lattice angs not supported 168 if '$lattice' in line: 169 # 0.00000000000 5.17898186576 5.17898186576 170 for lattvec in lines[l+1:l+4]: 171 lattvecs.append([float(x) for x in lattvec.split()]) 172 l += 4 173 # Look for Cartesian coordinates. They can be in another file or embedded in control: 174 #1) $coord file=coord 175 #2) $coord 176 # 2.58949092075 2.58949092075 2.58949092075 si 177 elif '$coord' in line: 178 if line.strip() == '$coord': 179 # Embdedded coordinates. 180 l += 1 181 while l < len(lines): 182 atom = lines[l].split() 183 if len(atom) == 4: 184 coords.append([float(x) for x in atom[0:3]]) 185 aspecies.append(atom[3].title()) # Convert si to Si, c to C, etc. 186 l += 1 187 else: 188 # End of $coord, go back to the main while loop to interpret the current line 189 break 190 elif line.find('file=') > 6: 191 # Cross-reference to another file 192 coordfile = line.split('=')[1].strip() 193 f_coord = open(coordfile) 194 for coordline in f_coord: 195 # 2.58949092075 2.58949092075 2.58949092075 si 196 atom = coordline.split() 197 if len(atom) == 4: 198 coords.append([float(x) for x in atom[0:3]]) 199 aspecies.append(atom[3].title()) # Convert si to Si, c to C, etc. 200 f_coord.close() 201 l += 1 202 else: 203 # $coordinateupdate or invalid $coord line 204 l += 1 205 else: 206 l += 1 207 208 if len(lattvecs) == 3 and len(aspecies) > 0 and len(aspecies) == len(coords): 209 self._tags['lattice_vectors'] = lattvecs 210 self._tags['atomic_species'] = aspecies 211 self._tags['coordinates'] = coords 212 else: 213 print("TURBOMOLE-interface: Error parsing TURBOMOLE output file") 214 215 216if __name__ == '__main__': 217 from phonopy.structure.symmetry import Symmetry 218 cell = read_turbomole(sys.argv[1]) 219 symmetry = Symmetry(cell) 220 print("# %s" % symmetry.get_international_table()) 221