Coverage for /builds/ase/ase/ase/io/castep/__init__.py: 61.76%

1# fmt: off 

2 

3"""This module defines I/O routines with CASTEP files. 

4The key idea is that all function accept or return atoms objects. 

5CASTEP specific parameters will be returned through the <atoms>.calc 

6attribute. 

7""" 

8import os 

9import re 

10import warnings 

11from copy import deepcopy 

12from typing import List, Tuple 

13 

14import numpy as np 

15 

16import ase 

17 

18# independent unit management included here: 

19# When high accuracy is required, this allows to easily pin down 

20# unit conversion factors from different "unit definition systems" 

21# (CODATA1986 for ase-3.6.0.2515 vs CODATA2002 for CASTEP 5.01). 

22# 

23# ase.units in in ase-3.6.0.2515 is based on CODATA1986 

24import ase.units 

25from ase.constraints import FixAtoms, FixCartesian, FixedLine, FixedPlane 

26from ase.geometry.cell import cellpar_to_cell 

27from ase.io.castep.castep_reader import read_castep_castep 

28from ase.parallel import paropen 

29from ase.spacegroup import Spacegroup 

30from ase.utils import atoms_to_spglib_cell, reader, writer 

31 

32from .geom_md_ts import ( 

33 read_castep_geom, 

34 read_castep_md, 

35 write_castep_geom, 

36 write_castep_md, 

37) 

38 

39units_ase = { 

40 'hbar': ase.units._hbar * ase.units.J, 

41 'Eh': ase.units.Hartree, 

42 'kB': ase.units.kB, 

43 'a0': ase.units.Bohr, 

44 't0': ase.units._hbar * ase.units.J / ase.units.Hartree, 

45 'c': ase.units._c, 

46 'me': ase.units._me / ase.units._amu, 

47 'Pascal': 1.0 / ase.units.Pascal} 

48 

49# CODATA1986 (included herein for the sake of completeness) 

50# taken from 

51# http://physics.nist.gov/cuu/Archive/1986RMP.pdf 

52units_CODATA1986 = { 

53 'hbar': 6.5821220E-16, # eVs 

54 'Eh': 27.2113961, # eV 

55 'kB': 8.617385E-5, # eV/K 

56 'a0': 0.529177249, # A 

57 'c': 299792458, # m/s 

58 'e': 1.60217733E-19, # C 

59 'me': 5.485799110E-4} # u 

60 

61# CODATA2002: default in CASTEP 5.01 

62# (-> check in more recent CASTEP in case of numerical discrepancies?!) 

63# taken from 

64# http://physics.nist.gov/cuu/Document/all_2002.pdf 

65units_CODATA2002 = { 

66 'hbar': 6.58211915E-16, # eVs 

67 'Eh': 27.2113845, # eV 

68 'kB': 8.617343E-5, # eV/K 

69 'a0': 0.5291772108, # A 

70 'c': 299792458, # m/s 

71 'e': 1.60217653E-19, # C 

72 'me': 5.4857990945E-4} # u 

73 

74# (common) derived entries 

75for d in (units_CODATA1986, units_CODATA2002): 

76 d['t0'] = d['hbar'] / d['Eh'] # s 

77 d['Pascal'] = d['e'] * 1E30 # Pa 

78 

79 

80__all__ = [ 

81 # routines for the generic io function 

82 'read_castep_castep', 

83 'read_castep_cell', 

84 'read_castep_geom', 

85 'read_castep_md', 

86 'read_phonon', 

87 'read_castep_phonon', 

88 # additional reads that still need to be wrapped 

89 'read_param', 

90 'read_seed', 

91 # write that is already wrapped 

92 'write_castep_cell', 

93 'write_castep_geom', 

94 'write_castep_md', 

95 # param write - in principle only necessary in junction with the calculator 

96 'write_param'] 

97 

98 

99def write_freeform(fd, outputobj): 

100 """ 

101 Prints out to a given file a CastepInputFile or derived class, such as 

102 CastepCell or CastepParam. 

103 """ 

104 

105 options = outputobj._options 

106 

107 # Some keywords, if present, are printed in this order 

108 preferred_order = ['lattice_cart', 'lattice_abc', 

109 'positions_frac', 'positions_abs', 

110 'species_pot', 'symmetry_ops', # CELL file 

111 'task', 'cut_off_energy' # PARAM file 

112 ] 

113 

114 keys = outputobj.get_attr_dict().keys() 

115 # This sorts only the ones in preferred_order and leaves the rest 

116 # untouched 

117 keys = sorted(keys, key=lambda x: preferred_order.index(x) 

118 if x in preferred_order 

119 else len(preferred_order)) 

120 

121 for kw in keys: 

122 opt = options[kw] 

123 if opt.type.lower() == 'block': 

124 fd.write('%BLOCK {0}\n{1}\n%ENDBLOCK {0}\n\n'.format( 

125 kw.upper(), 

126 opt.value.strip('\n'))) 

127 else: 

128 fd.write(f'{kw.upper()}: {opt.value}\n') 

129 

130 

131@writer 

132def write_castep_cell(fd, atoms, positions_frac=False, force_write=False, 

133 precision=6, magnetic_moments=None, 

134 castep_cell=None): 

135 """ 

136 This CASTEP export function write minimal information to 

137 a .cell file. If the atoms object is a trajectory, it will 

138 take the last image. 

139 

140 Note that function has been altered in order to require a filedescriptor 

141 rather than a filename. This allows to use the more generic write() 

142 function from formats.py 

143 

144 Note that the "force_write" keywords has no effect currently. 

145 

146 Arguments: 

147 

148 positions_frac: boolean. If true, positions are printed as fractional 

149 rather than absolute. Default is false. 

150 castep_cell: if provided, overrides the existing CastepCell object in 

151 the Atoms calculator 

152 precision: number of digits to which lattice and positions are printed 

153 magnetic_moments: if None, no SPIN values are initialised. 

154 If 'initial', the values from 

155 get_initial_magnetic_moments() are used. 

156 If 'calculated', the values from 

157 get_magnetic_moments() are used. 

158 If an array of the same length as the atoms object, 

159 its contents will be used as magnetic moments. 

160 """ 

161 

162 if isinstance(atoms, list): 

163 if len(atoms) > 1: 

164 atoms = atoms[-1] 

165 

166 # Header 

167 fd.write('# written by ASE\n\n') 

168 

169 # To write this we simply use the existing Castep calculator, or create 

170 # one 

171 from ase.calculators.castep import Castep, CastepCell 

172 

173 try: 

174 has_cell = isinstance(atoms.calc.cell, CastepCell) 

175 except AttributeError: 

176 has_cell = False 

177 

178 if has_cell: 

179 cell = deepcopy(atoms.calc.cell) 

180 else: 

181 cell = Castep(keyword_tolerance=2).cell 

182 

183 # Write lattice 

184 fformat = f'%{precision + 3}.{precision}f' 

185 cell_block_format = ' '.join([fformat] * 3) 

186 cell.lattice_cart = [cell_block_format % tuple(line) 

187 for line in atoms.get_cell()] 

188 

189 if positions_frac: 

190 pos_keyword = 'positions_frac' 

191 positions = atoms.get_scaled_positions() 

192 else: 

193 pos_keyword = 'positions_abs' 

194 positions = atoms.get_positions() 

195 

196 if atoms.has('castep_custom_species'): 

197 elems = atoms.get_array('castep_custom_species') 

198 else: 

199 elems = atoms.get_chemical_symbols() 

200 if atoms.has('masses'): 

201 

202 from ase.data import atomic_masses 

203 masses = atoms.get_array('masses') 

204 custom_masses = {} 

205 

206 for i, species in enumerate(elems): 

207 custom_mass = masses[i] 

208 

209 # build record of different masses for each species 

210 if species not in custom_masses: 

211 

212 # build dictionary of positions of all species with 

213 # same name and mass value ideally there should only 

214 # be one mass per species 

215 custom_masses[species] = {custom_mass: [i]} 

216 

217 # if multiple masses found for a species 

218 elif custom_mass not in custom_masses[species].keys(): 

219 

220 # if custom species were already manually defined raise an error 

221 if atoms.has('castep_custom_species'): 

222 raise ValueError( 

223 "Could not write custom mass block for {0}. \n" 

224 "Custom mass was set ({1}), but an inconsistent set of " 

225 "castep_custom_species already defines " 

226 "({2}) for {0}. \n" 

227 "If using both features, ensure that " 

228 "each species type in " 

229 "atoms.arrays['castep_custom_species'] " 

230 "has consistent mass values and that each atom " 

231 "with non-standard " 

232 "mass belongs to a custom species type." 

233 "".format( 

234 species, custom_mass, list( 

235 custom_masses[species].keys())[0])) 

236 

237 # append mass to create custom species later 

238 else: 

239 custom_masses[species][custom_mass] = [i] 

240 else: 

241 custom_masses[species][custom_mass].append(i) 

242 

243 # create species_mass block 

244 mass_block = [] 

245 

246 for el, mass_dict in custom_masses.items(): 

247 

248 # ignore mass record that match defaults 

249 default = mass_dict.pop(atomic_masses[atoms.get_array( 

250 'numbers')[list(elems).index(el)]], None) 

251 if mass_dict: 

252 # no custom species need to be created 

253 if len(mass_dict) == 1 and not default: 

254 mass_block.append('{} {}'.format( 

255 el, list(mass_dict.keys())[0])) 

256 # for each custom mass, create new species and change names to 

257 # match in 'elems' list 

258 else: 

259 warnings.warn( 

260 'Custom mass specified for ' 

261 'standard species {}, creating custom species' 

262 .format(el)) 

263 

264 for i, vals in enumerate(mass_dict.items()): 

265 mass_val, idxs = vals 

266 custom_species_name = f"{el}:{i}" 

267 warnings.warn( 

268 'Creating custom species {} with mass {}'.format( 

269 custom_species_name, str(mass_dict))) 

270 for idx in idxs: 

271 elems[idx] = custom_species_name 

272 mass_block.append('{} {}'.format( 

273 custom_species_name, mass_val)) 

274 

275 cell.species_mass = mass_block 

276 

277 if atoms.has('castep_labels'): 

278 labels = atoms.get_array('castep_labels') 

279 else: 

280 labels = ['NULL'] * len(elems) 

281 

282 if str(magnetic_moments).lower() == 'initial': 

283 magmoms = atoms.get_initial_magnetic_moments() 

284 elif str(magnetic_moments).lower() == 'calculated': 

285 magmoms = atoms.get_magnetic_moments() 

286 elif np.array(magnetic_moments).shape == (len(elems),): 

287 magmoms = np.array(magnetic_moments) 

288 else: 

289 magmoms = [0] * len(elems) 

290 

291 pos_block = [] 

292 pos_block_format = '%s ' + cell_block_format 

293 

294 for i, el in enumerate(elems): 

295 xyz = positions[i] 

296 line = pos_block_format % tuple([el] + list(xyz)) 

297 # ADD other keywords if necessary 

298 if magmoms[i] != 0: 

299 line += f' SPIN={magmoms[i]} ' 

300 if labels[i].strip() not in ('NULL', ''): 

301 line += f' LABEL={labels[i]} ' 

302 pos_block.append(line) 

303 

304 setattr(cell, pos_keyword, pos_block) 

305 

306 constr_block = _make_block_ionic_constraints(atoms) 

307 if constr_block: 

308 cell.ionic_constraints = constr_block 

309 

310 write_freeform(fd, cell) 

311 

312 

313def _make_block_ionic_constraints(atoms: ase.Atoms) -> List[str]: 

314 constr_block: List[str] = [] 

315 species_indices = atoms.symbols.species_indices() 

316 for constr in atoms.constraints: 

317 if not _is_constraint_valid(constr, len(atoms)): 

318 continue 

319 for i in constr.index: 

320 symbol = atoms.get_chemical_symbols()[i] 

321 nis = species_indices[i] + 1 

322 if isinstance(constr, FixAtoms): 

323 for j in range(3): # constraint for all three directions 

324 ic = len(constr_block) + 1 

325 line = f'{ic:6d} {symbol:3s} {nis:3d} ' 

326 line += ['1 0 0', '0 1 0', '0 0 1'][j] 

327 constr_block.append(line) 

328 elif isinstance(constr, FixCartesian): 

329 for j, m in enumerate(constr.mask): 

330 if m == 0: # not constrained 

331 continue 

332 ic = len(constr_block) + 1 

333 line = f'{ic:6d} {symbol:3s} {nis:3d} ' 

334 line += ['1 0 0', '0 1 0', '0 0 1'][j] 

335 constr_block.append(line) 

336 elif isinstance(constr, FixedPlane): 

337 ic = len(constr_block) + 1 

338 line = f'{ic:6d} {symbol:3s} {nis:3d} ' 

339 line += ' '.join([str(d) for d in constr.dir]) 

340 constr_block.append(line) 

341 elif isinstance(constr, FixedLine): 

342 for direction in _calc_normal_vectors(constr): 

343 ic = len(constr_block) + 1 

344 line = f'{ic:6d} {symbol:3s} {nis:3d} ' 

345 line += ' '.join(str(_) for _ in direction) 

346 constr_block.append(line) 

347 return constr_block 

348 

349 

350def _is_constraint_valid(constraint, natoms: int) -> bool: 

351 supported_constraints = (FixAtoms, FixedPlane, FixedLine, FixCartesian) 

352 if not isinstance(constraint, supported_constraints): 

353 warnings.warn(f'{constraint} is not supported by ASE CASTEP, skipped') 

354 return False 

355 if any(i < 0 or i >= natoms for i in constraint.index): 

356 warnings.warn(f'{constraint} contains invalid indices, skipped') 

357 return False 

358 return True 

359 

360 

361def _calc_normal_vectors(constr: FixedLine) -> Tuple[np.ndarray, np.ndarray]: 

362 direction = constr.dir 

363 

364 i2, i1 = np.argsort(np.abs(direction))[1:] 

365 v1 = direction[i1] 

366 v2 = direction[i2] 

367 n1 = np.zeros(3) 

368 n1[i2] = v1 

369 n1[i1] = -v2 

370 n1 = n1 / np.linalg.norm(n1) 

371 

372 n2 = np.cross(direction, n1) 

373 n2 = n2 / np.linalg.norm(n2) 

374 

375 return n1, n2 

376 

377 

378def read_freeform(fd): 

379 """ 

380 Read a CASTEP freeform file (the basic format of .cell and .param files) 

381 and return keyword-value pairs as a dict (values are strings for single 

382 keywords and lists of strings for blocks). 

383 """ 

384 

385 from ase.io.castep.castep_input_file import CastepInputFile 

386 

387 inputobj = CastepInputFile(keyword_tolerance=2) 

388 

389 filelines = fd.readlines() 

390 

391 keyw = None 

392 read_block = False 

393 block_lines = None 

394 

395 for i, l in enumerate(filelines): 

396 

397 # Strip all comments, aka anything after a hash 

398 L = re.split(r'[#!;]', l, maxsplit=1)[0].strip() 

399 

400 if L == '': 

401 # Empty line... skip 

402 continue 

403 

404 lsplit = re.split(r'\s*[:=]*\s+', L, maxsplit=1) 

405 

406 if read_block: 

407 if lsplit[0].lower() == '%endblock': 

408 if len(lsplit) == 1 or lsplit[1].lower() != keyw: 

409 raise ValueError('Out of place end of block at ' 

410 'line %i in freeform file' % i + 1) 

411 else: 

412 read_block = False 

413 inputobj.__setattr__(keyw, block_lines) 

414 else: 

415 block_lines += [L] 

416 else: 

417 # Check the first word 

418 

419 # Is it a block? 

420 read_block = (lsplit[0].lower() == '%block') 

421 if read_block: 

422 if len(lsplit) == 1: 

423 raise ValueError(('Unrecognizable block at line %i ' 

424 'in io freeform file') % i + 1) 

425 else: 

426 keyw = lsplit[1].lower() 

427 else: 

428 keyw = lsplit[0].lower() 

429 

430 # Now save the value 

431 if read_block: 

432 block_lines = [] 

433 else: 

434 inputobj.__setattr__(keyw, ' '.join(lsplit[1:])) 

435 

436 return inputobj.get_attr_dict(types=True) 

437 

438 

439@reader 

440def read_castep_cell(fd, index=None, calculator_args={}, find_spg=False, 

441 units=units_CODATA2002): 

442 """Read a .cell file and return an atoms object. 

443 Any value found that does not fit the atoms API 

444 will be stored in the atoms.calc attribute. 

445 

446 By default, the Castep calculator will be tolerant and in the absence of a 

447 castep_keywords.json file it will just accept all keywords that aren't 

448 automatically parsed. 

449 """ 

450 

451 from ase.calculators.castep import Castep 

452 

453 cell_units = { # Units specifiers for CASTEP 

454 'bohr': units_CODATA2002['a0'], 

455 'ang': 1.0, 

456 'm': 1e10, 

457 'cm': 1e8, 

458 'nm': 10, 

459 'pm': 1e-2 

460 } 

461 

462 calc = Castep(**calculator_args) 

463 

464 if calc.cell.castep_version == 0 and calc._kw_tol < 3: 

465 # No valid castep_keywords.json was found 

466 warnings.warn( 

467 'read_cell: Warning - Was not able to validate CASTEP input. ' 

468 'This may be due to a non-existing ' 

469 '"castep_keywords.json" ' 

470 'file or a non-existing CASTEP installation. ' 

471 'Parsing will go on but keywords will not be ' 

472 'validated and may cause problems if incorrect during a CASTEP ' 

473 'run.') 

474 

475 celldict = read_freeform(fd) 

476 

477 def parse_blockunit(line_tokens, blockname): 

478 u = 1.0 

479 if len(line_tokens[0]) == 1: 

480 usymb = line_tokens[0][0].lower() 

481 u = cell_units.get(usymb, 1) 

482 if usymb not in cell_units: 

483 warnings.warn('read_cell: Warning - ignoring invalid ' 

484 'unit specifier in %BLOCK {} ' 

485 '(assuming Angstrom instead)'.format(blockname)) 

486 line_tokens = line_tokens[1:] 

487 return u, line_tokens 

488 

489 # Arguments to pass to the Atoms object at the end 

490 aargs = { 

491 'pbc': True 

492 } 

493 

494 # Start by looking for the lattice 

495 lat_keywords = [w in celldict for w in ('lattice_cart', 'lattice_abc')] 

496 if all(lat_keywords): 

497 warnings.warn('read_cell: Warning - two lattice blocks present in the' 

498 ' same file. LATTICE_ABC will be ignored') 

499 elif not any(lat_keywords): 

500 raise ValueError('Cell file must contain at least one between ' 

501 'LATTICE_ABC and LATTICE_CART') 

502 

503 if 'lattice_abc' in celldict: 

504 

505 lines = celldict.pop('lattice_abc')[0].split('\n') 

506 line_tokens = [line.split() for line in lines] 

507 

508 u, line_tokens = parse_blockunit(line_tokens, 'lattice_abc') 

509 

510 if len(line_tokens) != 2: 

511 warnings.warn('read_cell: Warning - ignoring additional ' 

512 'lines in invalid %BLOCK LATTICE_ABC') 

513 

514 abc = [float(p) * u for p in line_tokens[0][:3]] 

515 angles = [float(phi) for phi in line_tokens[1][:3]] 

516 

517 aargs['cell'] = cellpar_to_cell(abc + angles) 

518 

519 if 'lattice_cart' in celldict: 

520 

521 lines = celldict.pop('lattice_cart')[0].split('\n') 

522 line_tokens = [line.split() for line in lines] 

523 

524 u, line_tokens = parse_blockunit(line_tokens, 'lattice_cart') 

525 

526 if len(line_tokens) != 3: 

527 warnings.warn('read_cell: Warning - ignoring more than ' 

528 'three lattice vectors in invalid %BLOCK ' 

529 'LATTICE_CART') 

530 

531 aargs['cell'] = [[float(x) * u for x in lt[:3]] for lt in line_tokens] 

532 

533 # Now move on to the positions 

534 pos_keywords = [w in celldict 

535 for w in ('positions_abs', 'positions_frac')] 

536 

537 if all(pos_keywords): 

538 warnings.warn('read_cell: Warning - two lattice blocks present in the' 

539 ' same file. POSITIONS_FRAC will be ignored') 

540 del celldict['positions_frac'] 

541 elif not any(pos_keywords): 

542 raise ValueError('Cell file must contain at least one between ' 

543 'POSITIONS_FRAC and POSITIONS_ABS') 

544 

545 aargs['symbols'] = [] 

546 pos_type = 'positions' 

547 pos_block = celldict.pop('positions_abs', [None])[0] 

548 if pos_block is None: 

549 pos_type = 'scaled_positions' 

550 pos_block = celldict.pop('positions_frac', [None])[0] 

551 aargs[pos_type] = [] 

552 

553 lines = pos_block.split('\n') 

554 line_tokens = [line.split() for line in lines] 

555 

556 if 'scaled' not in pos_type: 

557 u, line_tokens = parse_blockunit(line_tokens, 'positions_abs') 

558 else: 

559 u = 1.0 

560 

561 # Here we extract all the possible additional info 

562 # These are marked by their type 

563 

564 add_info = { 

565 'SPIN': (float, 0.0), # (type, default) 

566 'MAGMOM': (float, 0.0), 

567 'LABEL': (str, 'NULL') 

568 } 

569 add_info_arrays = {k: [] for k in add_info} 

570 

571 def parse_info(raw_info): 

572 

573 re_keys = (r'({0})\s*[=:\s]{{1}}\s' 

574 r'*([^\s]*)').format('|'.join(add_info.keys())) 

575 # Capture all info groups 

576 info = re.findall(re_keys, raw_info) 

577 info = {g[0]: add_info[g[0]][0](g[1]) for g in info} 

578 return info 

579 

580 # Array for custom species (a CASTEP special thing) 

581 # Usually left unused 

582 custom_species = None 

583 

584 for tokens in line_tokens: 

585 # Now, process the whole 'species' thing 

586 spec_custom = tokens[0].split(':', 1) 

587 elem = spec_custom[0] 

588 if len(spec_custom) > 1 and custom_species is None: 

589 # Add it to the custom info! 

590 custom_species = list(aargs['symbols']) 

591 if custom_species is not None: 

592 custom_species.append(tokens[0]) 

593 aargs['symbols'].append(elem) 

594 aargs[pos_type].append([float(p) * u for p in tokens[1:4]]) 

595 # Now for the additional information 

596 info = ' '.join(tokens[4:]) 

597 info = parse_info(info) 

598 for k in add_info: 

599 add_info_arrays[k] += [info.get(k, add_info[k][1])] 

600 

601 # read in custom species mass 

602 if 'species_mass' in celldict: 

603 spec_list = custom_species if custom_species else aargs['symbols'] 

604 aargs['masses'] = [None for _ in spec_list] 

605 lines = celldict.pop('species_mass')[0].split('\n') 

606 line_tokens = [line.split() for line in lines] 

607 

608 if len(line_tokens[0]) == 1: 

609 if line_tokens[0][0].lower() not in ('amu', 'u'): 

610 raise ValueError( 

611 "unit specifier '{}' in %BLOCK SPECIES_MASS " 

612 "not recognised".format( 

613 line_tokens[0][0].lower())) 

614 line_tokens = line_tokens[1:] 

615 

616 for tokens in line_tokens: 

617 token_pos_list = [i for i, x in enumerate( 

618 spec_list) if x == tokens[0]] 

619 if len(token_pos_list) == 0: 

620 warnings.warn( 

621 'read_cell: Warning - ignoring unused ' 

622 'species mass {} in %BLOCK SPECIES_MASS'.format( 

623 tokens[0])) 

624 for idx in token_pos_list: 

625 aargs['masses'][idx] = tokens[1] 

626 

627 # Now on to the species potentials... 

628 if 'species_pot' in celldict: 

629 lines = celldict.pop('species_pot')[0].split('\n') 

630 line_tokens = [line.split() for line in lines] 

631 

632 for tokens in line_tokens: 

633 if len(tokens) == 1: 

634 # It's a library 

635 all_spec = (set(custom_species) if custom_species is not None 

636 else set(aargs['symbols'])) 

637 for s in all_spec: 

638 calc.cell.species_pot = (s, tokens[0]) 

639 else: 

640 calc.cell.species_pot = tuple(tokens[:2]) 

641 

642 # Ionic constraints 

643 raw_constraints = {} 

644 

645 if 'ionic_constraints' in celldict: 

646 lines = celldict.pop('ionic_constraints')[0].split('\n') 

647 line_tokens = [line.split() for line in lines] 

648 

649 for tokens in line_tokens: 

650 if len(tokens) != 6: 

651 continue 

652 _, species, nic, x, y, z = tokens 

653 # convert xyz to floats 

654 x = float(x) 

655 y = float(y) 

656 z = float(z) 

657 

658 nic = int(nic) 

659 if (species, nic) not in raw_constraints: 

660 raw_constraints[(species, nic)] = [] 

661 raw_constraints[(species, nic)].append(np.array( 

662 [x, y, z])) 

663 

664 # Symmetry operations 

665 if 'symmetry_ops' in celldict: 

666 lines = celldict.pop('symmetry_ops')[0].split('\n') 

667 line_tokens = [line.split() for line in lines] 

668 

669 # Read them in blocks of four 

670 blocks = np.array(line_tokens).astype(float) 

671 if (len(blocks.shape) != 2 or blocks.shape[1] != 3 

672 or blocks.shape[0] % 4 != 0): 

673 warnings.warn('Warning: could not parse SYMMETRY_OPS' 

674 ' block properly, skipping') 

675 else: 

676 blocks = blocks.reshape((-1, 4, 3)) 

677 rotations = blocks[:, :3] 

678 translations = blocks[:, 3] 

679 

680 # Regardless of whether we recognize them, store these 

681 calc.cell.symmetry_ops = (rotations, translations) 

682 

683 # Anything else that remains, just add it to the cell object: 

684 for k, (val, otype) in celldict.items(): 

685 try: 

686 if otype == 'block': 

687 val = val.split('\n') # Avoids a bug for one-line blocks 

688 calc.cell.__setattr__(k, val) 

689 except Exception as e: 

690 raise RuntimeError( 

691 f'Problem setting calc.cell.{k} = {val}: {e}') 

692 

693 # Get the relevant additional info 

694 aargs['magmoms'] = np.array(add_info_arrays['SPIN']) 

695 # SPIN or MAGMOM are alternative keywords 

696 aargs['magmoms'] = np.where(aargs['magmoms'] != 0, 

697 aargs['magmoms'], 

698 add_info_arrays['MAGMOM']) 

699 labels = np.array(add_info_arrays['LABEL']) 

700 

701 aargs['calculator'] = calc 

702 

703 atoms = ase.Atoms(**aargs) 

704 

705 # Spacegroup... 

706 if find_spg: 

707 # Try importing spglib 

708 try: 

709 import spglib 

710 except ImportError: 

711 warnings.warn('spglib not found installed on this system - ' 

712 'automatic spacegroup detection is not possible') 

713 spglib = None 

714 

715 if spglib is not None: 

716 symmd = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms)) 

717 atoms_spg = Spacegroup(int(symmd['number'])) 

718 atoms.info['spacegroup'] = atoms_spg 

719 

720 atoms.new_array('castep_labels', labels) 

721 if custom_species is not None: 

722 atoms.new_array('castep_custom_species', np.array(custom_species)) 

723 

724 fixed_atoms = [] 

725 constraints = [] 

726 index_dict = atoms.symbols.indices() 

727 for (species, nic), value in raw_constraints.items(): 

728 

729 absolute_nr = index_dict[species][nic - 1] 

730 if len(value) == 3: 

731 # Check if they are linearly independent 

732 if np.linalg.det(value) == 0: 

733 warnings.warn( 

734 'Error: Found linearly dependent constraints attached ' 

735 'to atoms %s' % 

736 (absolute_nr)) 

737 continue 

738 fixed_atoms.append(absolute_nr) 

739 elif len(value) == 2: 

740 direction = np.cross(value[0], value[1]) 

741 # Check if they are linearly independent 

742 if np.linalg.norm(direction) == 0: 

743 warnings.warn( 

744 'Error: Found linearly dependent constraints attached ' 

745 'to atoms %s' % 

746 (absolute_nr)) 

747 continue 

748 constraint = FixedLine(indices=absolute_nr, direction=direction) 

749 constraints.append(constraint) 

750 elif len(value) == 1: 

751 direction = np.array(value[0], dtype=float) 

752 constraint = FixedPlane(indices=absolute_nr, direction=direction) 

753 constraints.append(constraint) 

754 else: 

755 warnings.warn( 

756 f'Error: Found {len(value)} statements attached to atoms ' 

757 f'{absolute_nr}' 

758 ) 

759 

760 # we need to sort the fixed atoms list in order not to raise an assertion 

761 # error in FixAtoms 

762 if fixed_atoms: 

763 constraints.append(FixAtoms(indices=sorted(fixed_atoms))) 

764 if constraints: 

765 atoms.set_constraint(constraints) 

766 

767 atoms.calc.atoms = atoms 

768 atoms.calc.push_oldstate() 

769 

770 return atoms 

771 

772 

773def read_phonon(filename, index=None, read_vib_data=False, 

774 gamma_only=True, frequency_factor=None, 

775 units=units_CODATA2002): 

776 """ 

777 Wrapper function for the more generic read() functionality. 

778 

779 Note that this is function is intended to maintain backwards-compatibility 

780 only. For documentation see read_castep_phonon(). 

781 """ 

782 from ase.io import read 

783 

784 if read_vib_data: 

785 full_output = True 

786 else: 

787 full_output = False 

788 

789 return read(filename, index=index, format='castep-phonon', 

790 full_output=full_output, read_vib_data=read_vib_data, 

791 gamma_only=gamma_only, frequency_factor=frequency_factor, 

792 units=units) 

793 

794 

795def read_castep_phonon(fd, index=None, read_vib_data=False, 

796 gamma_only=True, frequency_factor=None, 

797 units=units_CODATA2002): 

798 """ 

799 Reads a .phonon file written by a CASTEP Phonon task and returns an atoms 

800 object, as well as the calculated vibrational data if requested. 

801 

802 Note that the index argument has no effect as of now. 

803 """ 

804 

805 # fd is closed by embracing read() routine 

806 lines = fd.readlines() 

807 

808 atoms = None 

809 cell = [] 

810 N = Nb = Nq = 0 

811 scaled_positions = [] 

812 symbols = [] 

813 masses = []