1# Copyright (C) 2014 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_scaled_positions_lines, check_forces, 38 get_drift_forces) 39from phonopy.units import Bohr 40from phonopy.cui.settings import fracval 41from phonopy.structure.atoms import PhonopyAtoms as Atoms 42from phonopy.structure.symmetry import Symmetry 43 44# Castep output contains blank lines surrounded by astericks which should be ignored. 45# This is not implemented in collect_forces() function from file_IO module. 46# Below is the reimplementation of the collect_forces() function for Castep with 47# only one new variable skipafterhook. Setting this variable to zero (default value) 48# is equivalent to original collect_forces() function. 49def collect_forces_castep(f, num_atom, hook, force_pos, word=None, skipafterhook=0): 50 for line in f: 51 if hook in line: 52 break 53 if (skipafterhook>0): 54 for i in range(skipafterhook): 55 f.readline() 56 57 forces = [] 58 for line in f: 59 if line.strip() == '': 60 continue 61 if word is not None: 62 if word not in line: 63 continue 64 65 elems = line.split() 66 if len(elems) > force_pos[2]: 67 try: 68 forces.append([float(elems[i]) for i in force_pos]) 69 except ValueError: 70 forces = [] 71 break 72 else: 73 return False 74 75 if len(forces) == num_atom: 76 break 77 78 return forces 79 80def parse_set_of_forces(num_atoms, forces_filenames, verbose=True): 81 hook = 'Cartesian components (eV/A)' 82 skipafterhook=3 83 is_parsed = True 84 force_sets = [] 85 for i, filename in enumerate(forces_filenames): 86 if verbose: 87 sys.stdout.write("%d. " % (i + 1)) 88 89 f = open(filename) 90 castep_forces = collect_forces_castep(f, num_atoms, hook, [3, 4, 5], skipafterhook=3) 91 92 if check_forces(castep_forces, num_atoms, filename, verbose=verbose): 93 drift_force = get_drift_forces(castep_forces, 94 filename=filename, 95 verbose=verbose) 96 force_sets.append(np.array(castep_forces) - drift_force) 97 else: 98 is_parsed = False 99 100 if is_parsed: 101 return force_sets 102 else: 103 return [] 104 105def read_castep(filename): 106 f_castep = open(filename) 107 castep_in = CastepIn(f_castep.readlines()) 108 f_castep.close() 109 tags = castep_in.get_tags() 110# 1st stage is to create Atoms object ignoring Spin polarization. General case. 111 cell = Atoms(cell=tags['lattice_vectors'], 112 symbols=tags['atomic_species'], 113 scaled_positions=tags['coordinates']) 114# Analyse spin states and add data to Atoms instance "cell" if ones exist 115 magmoms = tags['magnetic_moments'] 116 if magmoms is not None: 117 # Print out symmetry information for magnetic cases 118 # Original code from structure/symmetry.py 119 symmetry = Symmetry(cell, symprec=1e-5) 120 print("CASTEP-interface: Magnetic structure, number of operations without spin: %d" % 121 len(symmetry.get_symmetry_operations()['rotations'])) 122 print("CASTEP-interface: Spacegroup without spin: %s" % symmetry.get_international_table()) 123 124 cell.set_magnetic_moments(magmoms) 125 symmetry = Symmetry(cell, symprec=1e-5) 126 print("CASTEP-interface: Magnetic structure, number of operations with spin: %d" % 127 len(symmetry.get_symmetry_operations()['rotations'])) 128 print("") 129 130 return cell 131 132def write_castep(filename, cell): 133 with open(filename, 'w') as f: 134 f.write(get_castep_structure(cell)) 135 136def write_supercells_with_displacements(supercell, 137 cells_with_displacements, 138 ids, 139 pre_filename="supercell", 140 width=3): 141 write_castep("%s.cell" % pre_filename, supercell) 142 for i, cell in zip(ids, cells_with_displacements): 143 filename = "{pre_filename}-{0:0{width}}.cell".format( 144 i, pre_filename=pre_filename, width=width) 145 write_castep(filename, cell) 146 147 148def get_castep_structure(cell): 149 lines = "" 150 lines += '%BLOCK LATTICE_CART\n' 151 lines += ((" % 20.16f" * 3 + "\n") * 3) % tuple(cell.get_cell().ravel()) 152 lines += '%ENDBLOCK LATTICE_CART\n\n' 153 lines += '%BLOCK POSITIONS_FRAC\n' 154 magmoms = cell.get_magnetic_moments() 155 156 for i in range(len(cell.get_chemical_symbols())): 157 atpos="".join("% 12.10f " % ap for ap in cell.get_scaled_positions()[i]) 158# Spin polarized case 159 if (magmoms is not None) and (magmoms[i] != 0.0): 160 lines += "".join("%2s %s spin=% 5.2f\n" % (cell.get_chemical_symbols()[i], atpos, magmoms[i])) 161# No spin ordering 162 else: 163 lines += "".join("%2s %s\n" % (cell.get_chemical_symbols()[i], atpos)) 164 lines += '%ENDBLOCK POSITIONS_FRAC\n\n' 165 166 return lines 167 168class CastepIn: 169 def __init__(self, lines): 170 self._tags = {'lattice_vectors': None, 171 'atomic_species': None, 172 'coordinates': None, 173 'magnetic_moments': None} 174 175 self._values = None 176 self._collect(lines) 177 178 def get_tags(self): 179 return self._tags 180 181 def _collect(self, lines): 182 magmoms = [] 183 184 numspins = 0 185 units = 1.0 # Angstrom units 186 lattvecs = [] 187 aspecies = [] 188 coords = [] 189 isSpinPol=0 190# Read cell parameters 191 for line in lines: 192 if ('%BLOCK LATTICE_CART' in line.upper()): 193 indx=lines.index(line) 194 for i in range(indx+1,len(lines)): 195 if ('ENDBLOCK' in lines[i].upper()): 196 break 197 if ('BOHR' in lines[i].upper()): 198 # The lattice vector units is Bohr. Convertion needed. 199 units=0.529177211 200 elif (len(lines[i].split())>=3): 201 lattvecs.append([(float(lines[i].split()[j])*units) for j in range(3)]) 202 203# Read atomic positions 204 for line in lines: 205 if ('%BLOCK POSITIONS_FRAC' in line.upper()): 206 indx=lines.index(line) 207 for i in range(indx+1,len(lines)): 208 if ('ENDBLOCK' in lines[i].upper()): 209 break 210 if (len(lines[i].split())>=3): 211 aspecies.append(lines[i].split()[0]) 212 coords.append([float(lines[i].split()[j]) for j in (1,2,3)]) 213 # If there is magmetic spin 214 if (len(lines[i].split())>4) and ('SPIN' in lines[i].upper()): 215 magmoms.append(float(lines[i].upper().split('SPIN')[1].split('=')[1])) 216 isSpinPol=1 217 else: 218 magmoms.append(0.0) 219# Set magnetic tags if the structure with magnetic order 220 if (isSpinPol>0): 221 self._tags['magnetic_moments'] = magmoms 222 223 if len(lattvecs) == 3 and len(aspecies) > 0 and len(aspecies) == len(coords): 224 self._tags['lattice_vectors'] = lattvecs 225 self._tags['atomic_species'] = aspecies 226 self._tags['coordinates'] = coords 227 else: 228 print("CASTEP-interface: Error parsing CASTEP .cell file") 229 230# ---------------------------------------------------------------------- 231