Coverage for 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 ( 

31 atoms_to_spglib_cell, 

32 reader, 

33 spglib_new_errorhandling, 

34 writer, 

35) 

36 

37from .geom_md_ts import ( 

38 read_castep_geom, 

39 read_castep_md, 

40 write_castep_geom, 

41 write_castep_md, 

42) 

43 

44units_ase = { 

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

46 'Eh': ase.units.Hartree, 

47 'kB': ase.units.kB, 

48 'a0': ase.units.Bohr, 

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

50 'c': ase.units._c, 

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

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

53 

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

55# taken from 

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

57units_CODATA1986 = { 

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

59 'Eh': 27.2113961, # eV 

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

61 'a0': 0.529177249, # A 

62 'c': 299792458, # m/s 

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

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

65 

66# CODATA2002: default in CASTEP 5.01 

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

68# taken from 

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

70units_CODATA2002 = { 

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

72 'Eh': 27.2113845, # eV 

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

74 'a0': 0.5291772108, # A 

75 'c': 299792458, # m/s 

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

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

78 

79# (common) derived entries 

80for d in (units_CODATA1986, units_CODATA2002): 

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

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

83 

84 

85__all__ = [ 

86 # routines for the generic io function 

87 'read_castep_castep', 

88 'read_castep_cell', 

89 'read_castep_geom', 

90 'read_castep_md', 

91 'read_phonon', 

92 'read_castep_phonon', 

93 # additional reads that still need to be wrapped 

94 'read_param', 

95 'read_seed', 

96 # write that is already wrapped 

97 'write_castep_cell', 

98 'write_castep_geom', 

99 'write_castep_md', 

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

101 'write_param'] 

102 

103 

104def write_freeform(fd, outputobj): 

105 """ 

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

107 CastepCell or CastepParam. 

108 """ 

109 

110 options = outputobj._options 

111 

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

113 preferred_order = ['lattice_cart', 'lattice_abc', 

114 'positions_frac', 'positions_abs', 

115 'species_pot', 'symmetry_ops', # CELL file 

116 'task', 'cut_off_energy' # PARAM file 

117 ] 

118 

119 keys = outputobj.get_attr_dict().keys() 

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

121 # untouched 

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

123 if x in preferred_order 

124 else len(preferred_order)) 

125 

126 for kw in keys: 

127 opt = options[kw] 

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

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

130 kw.upper(), 

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

132 else: 

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

134 

135 

136@writer 

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

138 precision=6, magnetic_moments=None, 

139 castep_cell=None): 

140 """ 

141 This CASTEP export function write minimal information to 

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

143 take the last image. 

144 

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

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

147 function from formats.py 

148 

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

150 

151 Arguments: 

152 

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

154 rather than absolute. Default is false. 

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

156 the Atoms calculator 

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

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

159 If 'initial', the values from 

160 get_initial_magnetic_moments() are used. 

161 If 'calculated', the values from 

162 get_magnetic_moments() are used. 

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

164 its contents will be used as magnetic moments. 

165 """ 

166 

167 if isinstance(atoms, list): 

168 if len(atoms) > 1: 

169 atoms = atoms[-1] 

170 

171 # Header 

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

173 

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

175 # one 

176 from ase.calculators.castep import Castep, CastepCell 

177 

178 try: 

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

180 except AttributeError: 

181 has_cell = False 

182 

183 if has_cell: 

184 cell = deepcopy(atoms.calc.cell) 

185 else: 

186 cell = Castep(keyword_tolerance=2).cell 

187 

188 # Write lattice 

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

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

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

192 for line in atoms.get_cell()] 

193 

194 if positions_frac: 

195 pos_keyword = 'positions_frac' 

196 positions = atoms.get_scaled_positions() 

197 else: 

198 pos_keyword = 'positions_abs' 

199 positions = atoms.get_positions() 

200 

201 if atoms.has('castep_custom_species'): 

202 elems = atoms.get_array('castep_custom_species') 

203 else: 

204 elems = atoms.get_chemical_symbols() 

205 if atoms.has('masses'): 

206 

207 from ase.data import atomic_masses 

208 masses = atoms.get_array('masses') 

209 custom_masses = {} 

210 

211 for i, species in enumerate(elems): 

212 custom_mass = masses[i] 

213 

214 # build record of different masses for each species 

215 if species not in custom_masses: 

216 

217 # build dictionary of positions of all species with 

218 # same name and mass value ideally there should only 

219 # be one mass per species 

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

221 

222 # if multiple masses found for a species 

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

224 

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

226 if atoms.has('castep_custom_species'): 

227 raise ValueError( 

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

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

230 "castep_custom_species already defines " 

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

232 "If using both features, ensure that " 

233 "each species type in " 

234 "atoms.arrays['castep_custom_species'] " 

235 "has consistent mass values and that each atom " 

236 "with non-standard " 

237 "mass belongs to a custom species type." 

238 "".format( 

239 species, custom_mass, list( 

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

241 

242 # append mass to create custom species later 

243 else: 

244 custom_masses[species][custom_mass] = [i] 

245 else: 

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

247 

248 # create species_mass block 

249 mass_block = [] 

250 

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

252 

253 # ignore mass record that match defaults 

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

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

256 if mass_dict: 

257 # no custom species need to be created 

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

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

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

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

262 # match in 'elems' list 

263 else: 

264 warnings.warn( 

265 'Custom mass specified for ' 

266 'standard species {}, creating custom species' 

267 .format(el)) 

268 

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

270 mass_val, idxs = vals 

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

272 warnings.warn( 

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

274 custom_species_name, str(mass_dict))) 

275 for idx in idxs: 

276 elems[idx] = custom_species_name 

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

278 custom_species_name, mass_val)) 

279 

280 cell.species_mass = mass_block 

281 

282 if atoms.has('castep_labels'): 

283 labels = atoms.get_array('castep_labels') 

284 else: 

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

286 

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

288 magmoms = atoms.get_initial_magnetic_moments() 

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

290 magmoms = atoms.get_magnetic_moments() 

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

292 magmoms = np.array(magnetic_moments) 

293 else: 

294 magmoms = [0] * len(elems) 

295 

296 pos_block = [] 

297 pos_block_format = '%s ' + cell_block_format 

298 

299 for i, el in enumerate(elems): 

300 xyz = positions[i] 

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

302 # ADD other keywords if necessary 

303 if magmoms[i] != 0: 

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

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

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

307 pos_block.append(line) 

308 

309 setattr(cell, pos_keyword, pos_block) 

310 

311 constr_block = _make_block_ionic_constraints(atoms) 

312 if constr_block: 

313 cell.ionic_constraints = constr_block 

314 

315 write_freeform(fd, cell) 

316 

317 

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

319 constr_block: List[str] = [] 

320 species_indices = atoms.symbols.species_indices() 

321 for constr in atoms.constraints: 

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

323 continue 

324 for i in constr.index: 

325 symbol = atoms.get_chemical_symbols()[i] 

326 nis = species_indices[i] + 1 

327 if isinstance(constr, FixAtoms): 

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

329 ic = len(constr_block) + 1 

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

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

332 constr_block.append(line) 

333 elif isinstance(constr, FixCartesian): 

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

335 if m == 0: # not constrained 

336 continue 

337 ic = len(constr_block) + 1 

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

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

340 constr_block.append(line) 

341 elif isinstance(constr, FixedPlane): 

342 ic = len(constr_block) + 1 

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

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

345 constr_block.append(line) 

346 elif isinstance(constr, FixedLine): 

347 for direction in _calc_normal_vectors(constr): 

348 ic = len(constr_block) + 1 

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

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

351 constr_block.append(line) 

352 return constr_block 

353 

354 

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

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

357 if not isinstance(constraint, supported_constraints): 

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

359 return False 

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

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

362 return False 

363 return True 

364 

365 

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

367 direction = constr.dir 

368 

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

370 v1 = direction[i1] 

371 v2 = direction[i2] 

372 n1 = np.zeros(3) 

373 n1[i2] = v1 

374 n1[i1] = -v2 

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

376 

377 n2 = np.cross(direction, n1) 

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

379 

380 return n1, n2 

381 

382 

383def read_freeform(fd): 

384 """ 

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

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

387 keywords and lists of strings for blocks). 

388 """ 

389 

390 from ase.io.castep.castep_input_file import CastepInputFile 

391 

392 inputobj = CastepInputFile(keyword_tolerance=2) 

393 

394 filelines = fd.readlines() 

395 

396 keyw = None 

397 read_block = False 

398 block_lines = None 

399 

400 for i, l in enumerate(filelines): 

401 

402 # Strip all comments, aka anything after a hash 

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

404 

405 if L == '': 

406 # Empty line... skip 

407 continue 

408 

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

410 

411 if read_block: 

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

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

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

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

416 else: 

417 read_block = False 

418 inputobj.__setattr__(keyw, block_lines) 

419 else: 

420 block_lines += [L] 

421 else: 

422 # Check the first word 

423 

424 # Is it a block? 

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

426 if read_block: 

427 if len(lsplit) == 1: 

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

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

430 else: 

431 keyw = lsplit[1].lower() 

432 else: 

433 keyw = lsplit[0].lower() 

434 

435 # Now save the value 

436 if read_block: 

437 block_lines = [] 

438 else: 

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

440 

441 return inputobj.get_attr_dict(types=True) 

442 

443 

444@reader 

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

446 units=units_CODATA2002): 

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

448 Any value found that does not fit the atoms API 

449 will be stored in the atoms.calc attribute. 

450 

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

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

453 automatically parsed. 

454 """ 

455 

456 from ase.calculators.castep import Castep 

457 

458 cell_units = { # Units specifiers for CASTEP 

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

460 'ang': 1.0, 

461 'm': 1e10, 

462 'cm': 1e8, 

463 'nm': 10, 

464 'pm': 1e-2 

465 } 

466 

467 calc = Castep(**calculator_args) 

468 

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

470 # No valid castep_keywords.json was found 

471 warnings.warn( 

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

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

474 '"castep_keywords.json" ' 

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

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

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

478 'run.') 

479 

480 celldict = read_freeform(fd) 

481 

482 def parse_blockunit(line_tokens, blockname): 

483 u = 1.0 

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

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

486 u = cell_units.get(usymb, 1) 

487 if usymb not in cell_units: 

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

489 'unit specifier in %BLOCK {} ' 

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

491 line_tokens = line_tokens[1:] 

492 return u, line_tokens 

493 

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

495 aargs = { 

496 'pbc': True 

497 } 

498 

499 # Start by looking for the lattice 

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

501 if all(lat_keywords): 

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

503 ' same file. LATTICE_ABC will be ignored') 

504 elif not any(lat_keywords): 

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

506 'LATTICE_ABC and LATTICE_CART') 

507 

508 if 'lattice_abc' in celldict: 

509 

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

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

512 

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

514 

515 if len(line_tokens) != 2: 

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

517 'lines in invalid %BLOCK LATTICE_ABC') 

518 

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

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

521 

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

523 

524 if 'lattice_cart' in celldict: 

525 

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

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

528 

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

530 

531 if len(line_tokens) != 3: 

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

533 'three lattice vectors in invalid %BLOCK ' 

534 'LATTICE_CART') 

535 

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

537 

538 # Now move on to the positions 

539 pos_keywords = [w in celldict 

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

541 

542 if all(pos_keywords): 

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

544 ' same file. POSITIONS_FRAC will be ignored') 

545 del celldict['positions_frac'] 

546 elif not any(pos_keywords): 

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

548 'POSITIONS_FRAC and POSITIONS_ABS') 

549 

550 aargs['symbols'] = [] 

551 pos_type = 'positions' 

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

553 if pos_block is None: 

554 pos_type = 'scaled_positions' 

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

556 aargs[pos_type] = [] 

557 

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

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

560 

561 if 'scaled' not in pos_type: 

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

563 else: 

564 u = 1.0 

565 

566 # Here we extract all the possible additional info 

567 # These are marked by their type 

568 

569 add_info = { 

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

571 'MAGMOM': (float, 0.0), 

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

573 } 

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

575 

576 def parse_info(raw_info): 

577 

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

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

580 # Capture all info groups 

581 info = re.findall(re_keys, raw_info) 

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

583 return info 

584 

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

586 # Usually left unused 

587 custom_species = None 

588 

589 for tokens in line_tokens: 

590 # Now, process the whole 'species' thing 

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

592 elem = spec_custom[0] 

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

594 # Add it to the custom info! 

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

596 if custom_species is not None: 

597 custom_species.append(tokens[0]) 

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

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

600 # Now for the additional information 

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

602 info = parse_info(info) 

603 for k in add_info: 

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

605 

606 # read in custom species mass 

607 if 'species_mass' in celldict: 

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

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

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

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

612 

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

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

615 raise ValueError( 

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

617 "not recognised".format( 

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

619 line_tokens = line_tokens[1:] 

620 

621 for tokens in line_tokens: 

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

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

624 if len(token_pos_list) == 0: 

625 warnings.warn( 

626 'read_cell: Warning - ignoring unused ' 

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

628 tokens[0])) 

629 for idx in token_pos_list: 

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

631 

632 # Now on to the species potentials... 

633 if 'species_pot' in celldict: 

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

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

636 

637 for tokens in line_tokens: 

638 if len(tokens) == 1: 

639 # It's a library 

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

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

642 for s in all_spec: 

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

644 else: 

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

646 

647 # Ionic constraints 

648 raw_constraints = {} 

649 

650 if 'ionic_constraints' in celldict: 

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

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

653 

654 for tokens in line_tokens: 

655 if len(tokens) != 6: 

656 continue 

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

658 # convert xyz to floats 

659 x = float(x) 

660 y = float(y) 

661 z = float(z) 

662 

663 nic = int(nic) 

664 if (species, nic) not in raw_constraints: 

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

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

667 [x, y, z])) 

668 

669 # Symmetry operations 

670 if 'symmetry_ops' in celldict: 

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

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

673 

674 # Read them in blocks of four 

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

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

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

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

679 ' block properly, skipping') 

680 else: 

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

682 rotations = blocks[:, :3] 

683 translations = blocks[:, 3] 

684 

685 # Regardless of whether we recognize them, store these 

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

687 

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

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

690 try: 

691 if otype == 'block': 

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

693 calc.cell.__setattr__(k, val) 

694 except Exception as e: 

695 raise RuntimeError( 

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

697 

698 # Get the relevant additional info 

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

700 # SPIN or MAGMOM are alternative keywords 

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

702 aargs['magmoms'], 

703 add_info_arrays['MAGMOM']) 

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

705 

706 aargs['calculator'] = calc 

707 

708 atoms = ase.Atoms(**aargs) 

709 

710 # Spacegroup... 

711 if find_spg: 

712 # Try importing spglib 

713 try: 

714 import spglib 

715 except ImportError: 

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

717 'automatic spacegroup detection is not possible') 

718 spglib = None 

719 

720 if spglib is not None: 

721 symmd = spglib_new_errorhandling(spglib.get_symmetry_dataset)( 

722 atoms_to_spglib_cell(atoms)) 

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

724 atoms.info['spacegroup'] = atoms_spg 

725 

726 atoms.new_array('castep_labels', labels) 

727 if custom_species is not None: 

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

729 

730 fixed_atoms = [] 

731 constraints = [] 

732 index_dict = atoms.symbols.indices() 

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

734 

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

736 if len(value) == 3: 

737 # Check if they are linearly independent 

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

739 warnings.warn( 

740 'Error: Found linearly dependent constraints attached ' 

741 'to atoms %s' % 

742 (absolute_nr)) 

743 continue 

744 fixed_atoms.append(absolute_nr) 

745 elif len(value) == 2: 

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

747 # Check if they are linearly independent 

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

749 warnings.warn( 

750 'Error: Found linearly dependent constraints attached ' 

751 'to atoms %s' % 

752 (absolute_nr)) 

753 continue 

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

755 constraints.append(constraint) 

756 elif len(value) == 1: 

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

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

759 constraints.append(constraint) 

760 else: 

761 warnings.warn( 

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

763 f'{absolute_nr}' 

764 ) 

765 

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

767 # error in FixAtoms 

768 if fixed_atoms: 

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

770 if constraints: 

771 atoms.set_constraint(constraints) 

772 

773 atoms.calc.atoms = atoms 

774 atoms.calc.push_oldstate() 

775 

776 return atoms 

777 

778 

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

780 gamma_only=True, frequency_factor=None, 

781 units=units_CODATA2002): 

782 """ 

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

784 

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

786 only. For documentation see read_castep_phonon(). 

787 """ 

788 from ase.io import read 

789 

790 if read_vib_data: 

791 full_output = True 

792 else: 

793 full_output = False 

794 

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

796 full_output=full_output, read_vib_data=read_vib_data, 

797 gamma_only=gamma_only, frequency_factor=frequency_factor, 

798 units=units) 

799 

800 

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

802 gamma_only=True, frequency_factor=None, 

803 units=units_CODATA2002): 

804 """ 

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

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

807 

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

809 """ 

810 

811 # fd is closed by embracing read() routine 

812 lines = fd.readlines() 

813 

814 atoms = None 

815 cell = []