Coverage for ase / io / castep / __init__.py: 74.15%

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 

12 

13import numpy as np 

14 

15import ase 

16 

17# independent unit management included here: 

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

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

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

21# 

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

23import ase.units 

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

25from ase.geometry.cell import cellpar_to_cell 

26from ase.io.castep.castep_reader import read_castep_castep 

27from ase.parallel import paropen 

28from ase.spacegroup import Spacegroup 

29from ase.utils import ( 

30 atoms_to_spglib_cell, 

31 reader, 

32 spglib_new_errorhandling, 

33 writer, 

34) 

35 

36from .geom_md_ts import ( 

37 read_castep_geom, 

38 read_castep_md, 

39 write_castep_geom, 

40 write_castep_md, 

41) 

42 

43units_ase = { 

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

45 'Eh': ase.units.Hartree, 

46 'kB': ase.units.kB, 

47 'a0': ase.units.Bohr, 

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

49 'c': ase.units._c, 

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

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

52 

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

54# taken from 

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

56units_CODATA1986 = { 

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

58 'Eh': 27.2113961, # eV 

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

60 'a0': 0.529177249, # A 

61 'c': 299792458, # m/s 

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

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

64 

65# CODATA2002: default in CASTEP 5.01 

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

67# taken from 

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

69units_CODATA2002 = { 

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

71 'Eh': 27.2113845, # eV 

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

73 'a0': 0.5291772108, # A 

74 'c': 299792458, # m/s 

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

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

77 

78# (common) derived entries 

79for d in (units_CODATA1986, units_CODATA2002): 

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

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

82 

83 

84__all__ = [ 

85 # routines for the generic io function 

86 'read_castep_castep', 

87 'read_castep_cell', 

88 'read_castep_geom', 

89 'read_castep_md', 

90 'read_phonon', 

91 'read_castep_phonon', 

92 # additional reads that still need to be wrapped 

93 'read_param', 

94 'read_seed', 

95 # write that is already wrapped 

96 'write_castep_cell', 

97 'write_castep_geom', 

98 'write_castep_md', 

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

100 'write_param'] 

101 

102 

103def write_freeform(fd, outputobj): 

104 """ 

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

106 CastepCell or CastepParam. 

107 """ 

108 

109 options = outputobj._options 

110 

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

112 preferred_order = ['lattice_cart', 'lattice_abc', 

113 'positions_frac', 'positions_abs', 

114 'species_pot', 'symmetry_ops', # CELL file 

115 'task', 'cut_off_energy' # PARAM file 

116 ] 

117 

118 keys = outputobj.get_attr_dict().keys() 

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

120 # untouched 

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

122 if x in preferred_order 

123 else len(preferred_order)) 

124 

125 for kw in keys: 

126 opt = options[kw] 

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

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

129 kw.upper(), 

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

131 else: 

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

133 

134 

135@writer 

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

137 precision=6, magnetic_moments=None, 

138 castep_cell=None): 

139 """ 

140 This CASTEP export function write minimal information to 

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

142 take the last image. 

143 

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

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

146 function from formats.py 

147 

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

149 

150 Arguments: 

151 

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

153 rather than absolute. Default is false. 

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

155 the Atoms calculator 

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

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

158 If 'initial', the values from 

159 get_initial_magnetic_moments() are used. 

160 If 'calculated', the values from 

161 get_magnetic_moments() are used. 

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

163 its contents will be used as magnetic moments. 

164 """ 

165 

166 if isinstance(atoms, list): 

167 if len(atoms) > 1: 

168 atoms = atoms[-1] 

169 

170 # Header 

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

172 

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

174 # one 

175 from ase.calculators.castep import Castep, CastepCell 

176 

177 try: 

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

179 except AttributeError: 

180 has_cell = False 

181 

182 if has_cell: 

183 cell = deepcopy(atoms.calc.cell) 

184 else: 

185 cell = Castep(keyword_tolerance=2).cell 

186 

187 # Write lattice 

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

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

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

191 for line in atoms.get_cell()] 

192 

193 if positions_frac: 

194 pos_keyword = 'positions_frac' 

195 positions = atoms.get_scaled_positions() 

196 else: 

197 pos_keyword = 'positions_abs' 

198 positions = atoms.get_positions() 

199 

200 if atoms.has('castep_custom_species'): 

201 elems = atoms.get_array('castep_custom_species') 

202 else: 

203 elems = atoms.get_chemical_symbols() 

204 if atoms.has('masses'): 

205 

206 from ase.data import atomic_masses 

207 masses = atoms.get_array('masses') 

208 custom_masses = {} 

209 

210 for i, species in enumerate(elems): 

211 custom_mass = masses[i] 

212 

213 # build record of different masses for each species 

214 if species not in custom_masses: 

215 

216 # build dictionary of positions of all species with 

217 # same name and mass value ideally there should only 

218 # be one mass per species 

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

220 

221 # if multiple masses found for a species 

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

223 

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

225 if atoms.has('castep_custom_species'): 

226 raise ValueError( 

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

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

229 "castep_custom_species already defines " 

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

231 "If using both features, ensure that " 

232 "each species type in " 

233 "atoms.arrays['castep_custom_species'] " 

234 "has consistent mass values and that each atom " 

235 "with non-standard " 

236 "mass belongs to a custom species type." 

237 "".format( 

238 species, custom_mass, list( 

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

240 

241 # append mass to create custom species later 

242 else: 

243 custom_masses[species][custom_mass] = [i] 

244 else: 

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

246 

247 # create species_mass block 

248 mass_block = [] 

249 

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

251 

252 # ignore mass record that match defaults 

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

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

255 if mass_dict: 

256 # no custom species need to be created 

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

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

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

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

261 # match in 'elems' list 

262 else: 

263 warnings.warn( 

264 'Custom mass specified for ' 

265 'standard species {}, creating custom species' 

266 .format(el)) 

267 

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

269 mass_val, idxs = vals 

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

271 warnings.warn( 

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

273 custom_species_name, str(mass_dict))) 

274 for idx in idxs: 

275 elems[idx] = custom_species_name 

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

277 custom_species_name, mass_val)) 

278 

279 cell.species_mass = mass_block 

280 

281 if atoms.has('castep_labels'): 

282 labels = atoms.get_array('castep_labels') 

283 else: 

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

285 

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

287 magmoms = atoms.get_initial_magnetic_moments() 

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

289 magmoms = atoms.get_magnetic_moments() 

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

291 magmoms = np.array(magnetic_moments) 

292 else: 

293 magmoms = [0] * len(elems) 

294 

295 pos_block = [] 

296 pos_block_format = '%s ' + cell_block_format 

297 

298 for i, el in enumerate(elems): 

299 xyz = positions[i] 

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

301 # ADD other keywords if necessary 

302 if magmoms[i] != 0: 

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

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

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

306 pos_block.append(line) 

307 

308 setattr(cell, pos_keyword, pos_block) 

309 

310 constr_block = _make_block_ionic_constraints(atoms) 

311 if constr_block: 

312 cell.ionic_constraints = constr_block 

313 

314 write_freeform(fd, cell) 

315 

316 

317def _make_block_ionic_constraints(atoms: ase.Atoms) -> list[str]: 

318 constr_block: list[str] = [] 

319 species_indices = atoms.symbols.species_indices() 

320 for constr in atoms.constraints: 

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

322 continue 

323 for i in constr.index: 

324 symbol = atoms.get_chemical_symbols()[i] 

325 nis = species_indices[i] + 1 

326 if isinstance(constr, FixAtoms): 

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

328 ic = len(constr_block) + 1 

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

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

331 constr_block.append(line) 

332 elif isinstance(constr, FixCartesian): 

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

334 if m == 0: # not constrained 

335 continue 

336 ic = len(constr_block) + 1 

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

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

339 constr_block.append(line) 

340 elif isinstance(constr, FixedPlane): 

341 ic = len(constr_block) + 1 

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

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

344 constr_block.append(line) 

345 elif isinstance(constr, FixedLine): 

346 for direction in _calc_normal_vectors(constr): 

347 ic = len(constr_block) + 1 

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

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

350 constr_block.append(line) 

351 return constr_block 

352 

353 

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

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

356 if not isinstance(constraint, supported_constraints): 

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

358 return False 

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

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

361 return False 

362 return True 

363 

364 

365def _calc_normal_vectors(constr: FixedLine) -> tuple[np.ndarray, np.ndarray]: 

366 direction = constr.dir 

367 

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

369 v1 = direction[i1] 

370 v2 = direction[i2] 

371 n1 = np.zeros(3) 

372 n1[i2] = v1 

373 n1[i1] = -v2 

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

375 

376 n2 = np.cross(direction, n1) 

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

378 

379 return n1, n2 

380 

381 

382def read_freeform(fd): 

383 """ 

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

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

386 keywords and lists of strings for blocks). 

387 """ 

388 

389 from ase.io.castep.castep_input_file import CastepInputFile 

390 

391 inputobj = CastepInputFile(keyword_tolerance=2) 

392 

393 filelines = fd.readlines() 

394 

395 keyw = None 

396 read_block = False 

397 block_lines = None 

398 

399 for i, l in enumerate(filelines): 

400 

401 # Strip all comments, aka anything after a hash 

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

403 

404 if L == '': 

405 # Empty line... skip 

406 continue 

407 

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

409 

410 if read_block: 

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

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

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

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

415 else: 

416 read_block = False 

417 inputobj.__setattr__(keyw, block_lines) 

418 else: 

419 block_lines += [L] 

420 else: 

421 # Check the first word 

422 

423 # Is it a block? 

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

425 if read_block: 

426 if len(lsplit) == 1: 

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

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

429 else: 

430 keyw = lsplit[1].lower() 

431 else: 

432 keyw = lsplit[0].lower() 

433 

434 # Now save the value 

435 if read_block: 

436 block_lines = [] 

437 else: 

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

439 

440 return inputobj.get_attr_dict(types=True) 

441 

442 

443@reader 

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

445 units=units_CODATA2002): 

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

447 Any value found that does not fit the atoms API 

448 will be stored in the atoms.calc attribute. 

449 

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

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

452 automatically parsed. 

453 """ 

454 

455 from ase.calculators.castep import Castep 

456 

457 cell_units = { # Units specifiers for CASTEP 

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

459 'ang': 1.0, 

460 'm': 1e10, 

461 'cm': 1e8, 

462 'nm': 10, 

463 'pm': 1e-2 

464 } 

465 

466 calc = Castep(**calculator_args) 

467 

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

469 # No valid castep_keywords.json was found 

470 warnings.warn( 

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

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

473 '"castep_keywords.json" ' 

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

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

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

477 'run.') 

478 

479 celldict = read_freeform(fd) 

480 

481 def parse_blockunit(line_tokens, blockname): 

482 u = 1.0 

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

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

485 u = cell_units.get(usymb, 1) 

486 if usymb not in cell_units: 

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

488 'unit specifier in %BLOCK {} ' 

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

490 line_tokens = line_tokens[1:] 

491 return u, line_tokens 

492 

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

494 aargs = { 

495 'pbc': True 

496 } 

497 

498 # Start by looking for the lattice 

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

500 if all(lat_keywords): 

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

502 ' same file. LATTICE_ABC will be ignored') 

503 elif not any(lat_keywords): 

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

505 'LATTICE_ABC and LATTICE_CART') 

506 

507 if 'lattice_abc' in celldict: 

508 

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

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

511 

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

513 

514 if len(line_tokens) != 2: 

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

516 'lines in invalid %BLOCK LATTICE_ABC') 

517 

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

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

520 

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

522 

523 if 'lattice_cart' in celldict: 

524 

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

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

527 

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

529 

530 if len(line_tokens) != 3: 

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

532 'three lattice vectors in invalid %BLOCK ' 

533 'LATTICE_CART') 

534 

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

536 

537 # Now move on to the positions 

538 pos_keywords = [w in celldict 

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

540 

541 if all(pos_keywords): 

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

543 ' same file. POSITIONS_FRAC will be ignored') 

544 del celldict['positions_frac'] 

545 elif not any(pos_keywords): 

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

547 'POSITIONS_FRAC and POSITIONS_ABS') 

548 

549 aargs['symbols'] = [] 

550 pos_type = 'positions' 

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

552 if pos_block is None: 

553 pos_type = 'scaled_positions' 

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

555 aargs[pos_type] = [] 

556 

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

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

559 

560 if 'scaled' not in pos_type: 

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

562 else: 

563 u = 1.0 

564 

565 # Here we extract all the possible additional info 

566 # These are marked by their type 

567 

568 add_info = { 

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

570 'MAGMOM': (float, 0.0), 

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

572 } 

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

574 

575 def parse_info(raw_info): 

576 

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

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

579 # Capture all info groups 

580 info = re.findall(re_keys, raw_info) 

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

582 return info 

583 

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

585 # Usually left unused 

586 custom_species = None 

587 

588 for tokens in line_tokens: 

589 # Now, process the whole 'species' thing 

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

591 elem = spec_custom[0] 

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

593 # Add it to the custom info! 

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

595 if custom_species is not None: 

596 custom_species.append(tokens[0]) 

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

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

599 # Now for the additional information 

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

601 info = parse_info(info) 

602 for k in add_info: 

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

604 

605 # read in custom species mass 

606 if 'species_mass' in celldict: 

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

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

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

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

611 

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

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

614 raise ValueError( 

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

616 "not recognised".format( 

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

618 line_tokens = line_tokens[1:] 

619 

620 for tokens in line_tokens: 

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

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

623 if len(token_pos_list) == 0: 

624 warnings.warn( 

625 'read_cell: Warning - ignoring unused ' 

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

627 tokens[0])) 

628 for idx in token_pos_list: 

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

630 

631 # Now on to the species potentials... 

632 if 'species_pot' in celldict: 

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

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

635 

636 for tokens in line_tokens: 

637 if len(tokens) == 1: 

638 # It's a library 

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

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

641 for s in all_spec: 

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

643 else: 

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

645 

646 # Ionic constraints 

647 raw_constraints = {} 

648 

649 if 'ionic_constraints' in celldict: 

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

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

652 

653 for tokens in line_tokens: 

654 if len(tokens) != 6: 

655 continue 

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

657 # convert xyz to floats 

658 x = float(x) 

659 y = float(y) 

660 z = float(z) 

661 

662 nic = int(nic) 

663 if (species, nic) not in raw_constraints: 

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

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

666 [x, y, z])) 

667 

668 # Symmetry operations 

669 if 'symmetry_ops' in celldict: 

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

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

672 

673 # Read them in blocks of four 

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

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

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

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

678 ' block properly, skipping') 

679 else: 

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

681 rotations = blocks[:, :3] 

682 translations = blocks[:, 3] 

683 

684 # Regardless of whether we recognize them, store these 

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

686 

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

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

689 try: 

690 if otype == 'block': 

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

692 calc.cell.__setattr__(k, val) 

693 except Exception as e: 

694 raise RuntimeError( 

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

696 

697 # Get the relevant additional info 

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

699 # SPIN or MAGMOM are alternative keywords 

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

701 aargs['magmoms'], 

702 add_info_arrays['MAGMOM']) 

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

704 

705 aargs['calculator'] = calc 

706 

707 atoms = ase.Atoms(**aargs) 

708 

709 # Spacegroup... 

710 if find_spg: 

711 # Try importing spglib 

712 try: 

713 import spglib 

714 except ImportError: 

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

716 'automatic spacegroup detection is not possible') 

717 spglib = None 

718 

719 if spglib is not None: 

720 symmd = spglib_new_errorhandling(spglib.get_symmetry_dataset)( 

721 atoms_to_spglib_cell(atoms)) 

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

723 atoms.info['spacegroup'] = atoms_spg 

724 

725 atoms.new_array('castep_labels', labels) 

726 if custom_species is not None: 

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

728 

729 fixed_atoms = [] 

730 constraints = [] 

731 index_dict = atoms.symbols.indices() 

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

733 

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

735 if len(value) == 3: 

736 # Check if they are linearly independent 

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

738 warnings.warn( 

739 'Error: Found linearly dependent constraints attached ' 

740 'to atoms %s' % 

741 (absolute_nr)) 

742 continue 

743 fixed_atoms.append(absolute_nr) 

744 elif len(value) == 2: 

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

746 # Check if they are linearly independent 

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

748 warnings.warn( 

749 'Error: Found linearly dependent constraints attached ' 

750 'to atoms %s' % 

751 (absolute_nr)) 

752 continue 

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

754 constraints.append(constraint) 

755 elif len(value) == 1: 

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

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

758 constraints.append(constraint) 

759 else: 

760 warnings.warn( 

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

762 f'{absolute_nr}' 

763 ) 

764 

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

766 # error in FixAtoms 

767 if fixed_atoms: 

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

769 if constraints: 

770 atoms.set_constraint(constraints) 

771 

772 atoms.calc.atoms = atoms 

773 atoms.calc.push_oldstate() 

774 

775 return atoms 

776 

777 

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

779 gamma_only=True, frequency_factor=None, 

780 units=units_CODATA2002): 

781 """ 

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

783 

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

785 only. For documentation see read_castep_phonon(). 

786 """ 

787 from ase.io import read 

788 

789 if read_vib_data: 

790 full_output = True 

791 else: 

792 full_output = False 

793 

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

795 full_output=full_output, read_vib_data=read_vib_data, 

796 gamma_only=gamma_only, frequency_factor=frequency_factor, 

797 units=units) 

798 

799 

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

801 gamma_only=True, frequency_factor=None, 

802 units=units_CODATA2002): 

803 """ 

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

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

806 

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

808 """ 

809 

810 # fd is closed by embracing read() routine 

811 lines = fd.readlines() 

812 

813 atoms = None 

814 cell = [] 

815 N = Nb = Nq = 0