Coverage for /builds/alexhroom/ase/ase/io/castep/__init_

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

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

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

4attribute.

5"""

6import os

7import re

8import warnings

9from copy import deepcopy

10from typing import List, Tuple

12import numpy as np

14import ase

15# independent unit management included here:

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

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

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

19#

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

21import ase.units

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

23from ase.geometry.cell import cellpar_to_cell

24from ase.io.castep.castep_reader import read_castep_castep

25from ase.parallel import paropen

26from ase.spacegroup import Spacegroup

27from ase.utils import atoms_to_spglib_cell, reader

29units_ase = {

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

31 'Eh': ase.units.Hartree,

32 'kB': ase.units.kB,

33 'a0': ase.units.Bohr,

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

35 'c': ase.units._c,

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

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

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

40# taken from

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

42units_CODATA1986 = {

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

44 'Eh': 27.2113961, # eV

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

46 'a0': 0.529177249, # A

47 'c': 299792458, # m/s

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

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

51# CODATA2002: default in CASTEP 5.01

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

53# taken from

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

55units_CODATA2002 = {

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

57 'Eh': 27.2113845, # eV

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

59 'a0': 0.5291772108, # A

60 'c': 299792458, # m/s

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

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

64# (common) derived entries

65for d in (units_CODATA1986, units_CODATA2002):

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

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

70__all__ = [

71 # routines for the generic io function

72 'read_castep_castep',

73 'read_cell',

74 'read_castep_cell',

75 'read_geom',

76 'read_castep_geom',

77 'read_phonon',

78 'read_castep_phonon',

79 # additional reads that still need to be wrapped

80 'read_md',

81 'read_param',

82 'read_seed',

83 # write that is already wrapped

84 'write_castep_cell',

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

86 'write_param']

89def write_freeform(fd, outputobj):

90 """

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

92 CastepCell or CastepParam.

93 """

95 options = outputobj._options

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

98 preferred_order = ['lattice_cart', 'lattice_abc',

99 'positions_frac', 'positions_abs',

100 'species_pot', 'symmetry_ops', # CELL file

101 'task', 'cut_off_energy' # PARAM file

102 ]

103

104 keys = outputobj.get_attr_dict().keys()

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

106 # untouched

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

108 if x in preferred_order

109 else len(preferred_order))

110

111 for kw in keys:

112 opt = options[kw]

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

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

115 kw.upper(),

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

117 else:

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

119

120

121def write_cell(filename, atoms, positions_frac=False, castep_cell=None,

122 force_write=False):

123 """

124 Wrapper function for the more generic write() functionality.

125

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

127 only.

128 """

129 from ase.io import write

130

131 write(filename, atoms, positions_frac=positions_frac,

132 castep_cell=castep_cell, force_write=force_write)

133

134

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

136 precision=6, magnetic_moments=None,

137 castep_cell=None):

138 """

139 This CASTEP export function write minimal information to

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

141 take the last image.

142

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

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

145 function from formats.py

146

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

148

149 Arguments:

150

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

152 rather than absolute. Default is false.

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

154 the Atoms calculator

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

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

157 If 'initial', the values from

158 get_initial_magnetic_moments() are used.

159 If 'calculated', the values from

160 get_magnetic_moments() are used.

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

162 its contents will be used as magnetic moments.

163 """

164

165 if isinstance(atoms, list):

166 if len(atoms) > 1:

167 atoms = atoms[-1]

168

169 # Header

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

171

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

173 # one

174 from ase.calculators.castep import Castep, CastepCell

175

176 try:

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

178 except AttributeError:

179 has_cell = False

180

181 if has_cell:

182 cell = deepcopy(atoms.calc.cell)

183 else:

184 cell = Castep(keyword_tolerance=2).cell

185

186 # Write lattice

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

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

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

190 for line in atoms.get_cell()]

191

192 if positions_frac:

193 pos_keyword = 'positions_frac'

194 positions = atoms.get_scaled_positions()

195 else:

196 pos_keyword = 'positions_abs'

197 positions = atoms.get_positions()

198

199 if atoms.has('castep_custom_species'):

200 elems = atoms.get_array('castep_custom_species')

201 else:

202 elems = atoms.get_chemical_symbols()

203 if atoms.has('masses'):

204

205 from ase.data import atomic_masses

206 masses = atoms.get_array('masses')

207 custom_masses = {}

208

209 for i, species in enumerate(elems):

210 custom_mass = masses[i]

211

212 # build record of different masses for each species

213 if species not in custom_masses:

214

215 # build dictionary of positions of all species with

216 # same name and mass value ideally there should only

217 # be one mass per species

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

219

220 # if multiple masses found for a species

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

222

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

224 if atoms.has('castep_custom_species'):

225 raise ValueError(

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

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

228 "castep_custom_species already defines "

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

230 "If using both features, ensure that "

231 "each species type in "

232 "atoms.arrays['castep_custom_species'] "

233 "has consistent mass values and that each atom "

234 "with non-standard "

235 "mass belongs to a custom species type."

236 "".format(

237 species, custom_mass, list(

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

239

240 # append mass to create custom species later

241 else:

242 custom_masses[species][custom_mass] = [i]

243 else:

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

245

246 # create species_mass block

247 mass_block = []

248

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

250

251 # ignore mass record that match defaults

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

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

254 if mass_dict:

255 # no custom species need to be created

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

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

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

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

260 # match in 'elems' list

261 else:

262 warnings.warn(

263 'Custom mass specified for '

264 'standard species {}, creating custom species'

265 .format(el))

266

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

268 mass_val, idxs = vals

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

270 warnings.warn(

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

272 custom_species_name, str(mass_dict)))

273 for idx in idxs:

274 elems[idx] = custom_species_name

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

276 custom_species_name, mass_val))

277

278 cell.species_mass = mass_block

279

280 if atoms.has('castep_labels'):

281 labels = atoms.get_array('castep_labels')

282 else:

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

284

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

286 magmoms = atoms.get_initial_magnetic_moments()

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

288 magmoms = atoms.get_magnetic_moments()

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

290 magmoms = np.array(magnetic_moments)

291 else:

292 magmoms = [0] * len(elems)

293

294 pos_block = []

295 pos_block_format = '%s ' + cell_block_format

296

297 for i, el in enumerate(elems):

298 xyz = positions[i]

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

300 # ADD other keywords if necessary

301 if magmoms[i] != 0:

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

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

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

305 pos_block.append(line)

306

307 setattr(cell, pos_keyword, pos_block)

308

309 constr_block = _make_block_ionic_constraints(atoms)

310 if constr_block:

311 cell.ionic_constraints = constr_block

312

313 write_freeform(fd, cell)

314

315

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

317 constr_block: List[str] = []

318 species_indices = atoms.symbols.species_indices()

319 for constr in atoms.constraints:

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

321 continue

322 for i in constr.index:

323 symbol = atoms.get_chemical_symbols()[i]

324 nis = species_indices[i] + 1

325 if isinstance(constr, FixAtoms):

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

327 ic = len(constr_block) + 1

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

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

330 constr_block.append(line)

331 elif isinstance(constr, FixCartesian):

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

333 if m == 0: # not constrained

334 continue

335 ic = len(constr_block) + 1

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

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

338 constr_block.append(line)

339 elif isinstance(constr, FixedPlane):

340 ic = len(constr_block) + 1

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

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

343 constr_block.append(line)

344 elif isinstance(constr, FixedLine):

345 for direction in _calc_normal_vectors(constr):

346 ic = len(constr_block) + 1

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

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

349 constr_block.append(line)

350 return constr_block

351

352

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

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

355 if not isinstance(constraint, supported_constraints):

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

357 return False

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

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

360 return False

361 return True

362

363

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

365 direction = constr.dir

366

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

368 v1 = direction[i1]

369 v2 = direction[i2]

370 n1 = np.zeros(3)

371 n1[i2] = v1

372 n1[i1] = -v2

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

374

375 n2 = np.cross(direction, n1)

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

377

378 return n1, n2

379

380

381def read_freeform(fd):

382 """

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

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

385 keywords and lists of strings for blocks).

386 """

387

388 from ase.io.castep.castep_input_file import CastepInputFile

389

390 inputobj = CastepInputFile(keyword_tolerance=2)

391

392 filelines = fd.readlines()

393

394 keyw = None

395 read_block = False

396 block_lines = None

397

398 for i, l in enumerate(filelines):

399

400 # Strip all comments, aka anything after a hash

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

402

403 if L == '':

404 # Empty line... skip

405 continue

406

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

408

409 if read_block:

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

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

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

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

414 else:

415 read_block = False

416 inputobj.__setattr__(keyw, block_lines)

417 else:

418 block_lines += [L]

419 else:

420 # Check the first word

421

422 # Is it a block?

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

424 if read_block:

425 if len(lsplit) == 1:

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

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

428 else:

429 keyw = lsplit[1].lower()

430 else:

431 keyw = lsplit[0].lower()

432

433 # Now save the value

434 if read_block:

435 block_lines = []

436 else:

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

438

439 return inputobj.get_attr_dict(types=True)

440

441

442def read_cell(filename, index=None):

443 """

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

445

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

447 only.

448 """

449 from ase.io import read

450 return read(filename, index=index, format='castep-cell')

451

452

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

454 units=units_CODATA2002):

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

456 Any value found that does not fit the atoms API

457 will be stored in the atoms.calc attribute.

458

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

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

461 automatically parsed.

462 """

463

464 from ase.calculators.castep import Castep

465

466 cell_units = { # Units specifiers for CASTEP

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

468 'ang': 1.0,

469 'm': 1e10,

470 'cm': 1e8,

471 'nm': 10,

472 'pm': 1e-2

473 }

474

475 calc = Castep(**calculator_args)

476

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

478 # No valid castep_keywords.json was found

479 warnings.warn(

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

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

482 '"castep_keywords.json" '

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

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

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

486 'run.')

487

488 celldict = read_freeform(fd)

489

490 def parse_blockunit(line_tokens, blockname):

491 u = 1.0

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

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

494 u = cell_units.get(usymb, 1)

495 if usymb not in cell_units:

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

497 'unit specifier in %BLOCK {} '

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

499 line_tokens = line_tokens[1:]

500 return u, line_tokens

501

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

503 aargs = {

504 'pbc': True

505 }

506

507 # Start by looking for the lattice

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

509 if all(lat_keywords):

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

511 ' same file. LATTICE_ABC will be ignored')

512 elif not any(lat_keywords):

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

514 'LATTICE_ABC and LATTICE_CART')

515

516 if 'lattice_abc' in celldict:

517

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

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

520

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

522

523 if len(line_tokens) != 2:

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

525 'lines in invalid %BLOCK LATTICE_ABC')

526

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

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

529

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

531

532 if 'lattice_cart' in celldict:

533

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

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

536

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

538

539 if len(line_tokens) != 3:

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

541 'three lattice vectors in invalid %BLOCK '

542 'LATTICE_CART')

543

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

545

546 # Now move on to the positions

547 pos_keywords = [w in celldict

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

549

550 if all(pos_keywords):

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

552 ' same file. POSITIONS_FRAC will be ignored')

553 del celldict['positions_frac']

554 elif not any(pos_keywords):

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

556 'POSITIONS_FRAC and POSITIONS_ABS')

557

558 aargs['symbols'] = []

559 pos_type = 'positions'

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

561 if pos_block is None:

562 pos_type = 'scaled_positions'

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

564 aargs[pos_type] = []

565

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

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

568

569 if 'scaled' not in pos_type:

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

571 else:

572 u = 1.0

573

574 # Here we extract all the possible additional info

575 # These are marked by their type

576

577 add_info = {

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

579 'MAGMOM': (float, 0.0),

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

581 }

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

583

584 def parse_info(raw_info):

585

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

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

588 # Capture all info groups

589 info = re.findall(re_keys, raw_info)

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

591 return info

592

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

594 # Usually left unused

595 custom_species = None

596

597 for tokens in line_tokens:

598 # Now, process the whole 'species' thing

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

600 elem = spec_custom[0]

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

602 # Add it to the custom info!

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

604 if custom_species is not None:

605 custom_species.append(tokens[0])

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

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

608 # Now for the additional information

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

610 info = parse_info(info)

611 for k in add_info:

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

613

614 # read in custom species mass

615 if 'species_mass' in celldict:

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

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

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

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

620

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

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

623 raise ValueError(

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

625 "not recognised".format(

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

627 line_tokens = line_tokens[1:]

628

629 for tokens in line_tokens:

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

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

632 if len(token_pos_list) == 0:

633 warnings.warn(

634 'read_cell: Warning - ignoring unused '

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

636 tokens[0]))

637 for idx in token_pos_list:

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

639

640 # Now on to the species potentials...

641 if 'species_pot' in celldict:

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

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

644

645 for tokens in line_tokens:

646 if len(tokens) == 1:

647 # It's a library

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

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

650 for s in all_spec:

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

652 else:

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

654

655 # Ionic constraints

656 raw_constraints = {}

657

658 if 'ionic_constraints' in celldict:

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

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

661

662 for tokens in line_tokens:

663 if len(tokens) != 6:

664 continue

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

666 # convert xyz to floats

667 x = float(x)

668 y = float(y)

669 z = float(z)

670

671 nic = int(nic)

672 if (species, nic) not in raw_constraints:

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

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

675 [x, y, z]))

676

677 # Symmetry operations

678 if 'symmetry_ops' in celldict:

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

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

681

682 # Read them in blocks of four

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

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

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

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

687 ' block properly, skipping')

688 else:

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

690 rotations = blocks[:, :3]

691 translations = blocks[:, 3]

692

693 # Regardless of whether we recognize them, store these

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

695

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

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

698 try:

699 if otype == 'block':

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

701 calc.cell.__setattr__(k, val)

702 except Exception as e:

703 raise RuntimeError(

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

705

706 # Get the relevant additional info

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

708 # SPIN or MAGMOM are alternative keywords

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

710 aargs['magmoms'],

711 add_info_arrays['MAGMOM'])

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

713

714 aargs['calculator'] = calc

715

716 atoms = ase.Atoms(**aargs)

717

718 # Spacegroup...

719 if find_spg:

720 # Try importing spglib

721 try:

722 import spglib

723 except ImportError:

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

725 'automatic spacegroup detection is not possible')

726 spglib = None

727

728 if spglib is not None:

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

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

731 atoms.info['spacegroup'] = atoms_spg

732

733 atoms.new_array('castep_labels', labels)

734 if custom_species is not None:

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

736

737 fixed_atoms = []

738 constraints = []

739 index_dict = atoms.symbols.indices()

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

741

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

743 if len(value) == 3:

744 # Check if they are linearly independent

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

746 warnings.warn(

747 'Error: Found linearly dependent constraints attached '

748 'to atoms %s' %

749 (absolute_nr))

750 continue

751 fixed_atoms.append(absolute_nr)

752 elif len(value) == 2:

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

754 # Check if they are linearly independent

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

756 warnings.warn(

757 'Error: Found linearly dependent constraints attached '

758 'to atoms %s' %

759 (absolute_nr))

760 continue

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

762 constraints.append(constraint)

763 elif len(value) == 1:

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

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

766 constraints.append(constraint)

767 else:

768 warnings.warn(

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

770 f'{absolute_nr}'

771 )

772

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

774 # error in FixAtoms

775 if fixed_atoms:

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

777 if constraints:

778 atoms.set_constraint(constraints)

779

780 atoms.calc.atoms = atoms

781 atoms.calc.push_oldstate()

782

783 return atoms

784

785

786def read_geom(filename, index=':', units=units_CODATA2002):

787 """

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

789

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

791 only. Keyword arguments will be passed to read_castep_geom().

792 """

793 from ase.io import read

794 return read(filename, index=index, format='castep-geom', units=units)

795

796

797def read_castep_geom(fd, index=None, units=units_CODATA2002):

798 """Reads a .geom file produced by the CASTEP GeometryOptimization task and

799 returns an atoms object.

800 The information about total free energy and forces of each atom for every

801 relaxation step will be stored for further analysis especially in a

802 single-point calculator.

803 Note that everything in the .geom file is in atomic units, which has

804 been conversed to commonly used unit angstrom(length) and eV (energy).

805

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

807

808 Contribution by Wei-Bing Zhang. Thanks!

809

810 Routine now accepts a filedescriptor in order to out-source the gz and

811 bz2 handling to formats.py. Note that there is a fallback routine

812 read_geom() that behaves like previous versions did.

813 """

814 from ase.calculators.singlepoint import SinglePointCalculator

815

816 # fd is closed by embracing read() routine

817 txt = fd.readlines()

818

819 traj = []

820

821 Hartree = units['Eh']

822 Bohr = units['a0']

823

824 # Yeah, we know that...

825 # print('N.B.: Energy in .geom file is not 0K extrapolated.')

826 for i, line in enumerate(txt):

827 if line.find('<-- E') > 0:

828 start_found = True

829 energy = float(line.split()[0]) * Hartree

830 cell = [x.split()[0:3] for x in txt[i + 1:i + 4]]

831 cell = np.array([[float(col) * Bohr for col in row] for row in

832 cell])

833 if line.find('<-- R') > 0 and start_found:

834 start_found = False

835 geom_start = i

836 for i, line in enumerate(txt[geom_start:]):

837 if line.find('<-- F') > 0:

838 geom_stop = i + geom_start

839 break

840 species = [line.split()[0] for line in

841 txt[geom_start:geom_stop]]

842 geom = np.array([[float(col) * Bohr for col in

843 line.split()[2:5]] for line in

844 txt[geom_start:geom_stop]])

845 forces = np.array([[float(col) * Hartree / Bohr for col in

846 line.split()[2:5]] for line in

847 txt[geom_stop:geom_stop

848 + (geom_stop - geom_start)]])

849 image = ase.Atoms(species, geom, cell=cell, pbc=True)

850 image.calc = SinglePointCalculator(

851 atoms=image, energy=energy, forces=forces)

852 traj.append(image)

853

854 if index is None:

855 return traj

856 else:

857 return traj[index]

858

859

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

861 gamma_only=True, frequency_factor=None,

862 units=units_CODATA2002):

863 """

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

865

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

867 only. For documentation see read_castep_phonon().

868 """

869 from ase.io import read

870

871 if read_vib_data:

872 full_output = True

873 else:

874 full_output = False

875

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

877 full_output=full_output, read_vib_data=read_vib_data,

878 gamma_only=gamma_only, frequency_factor=frequency_factor,

879 units=units)

880

881

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

883 gamma_only=True, frequency_factor=None,

884 units=units_CODATA2002):

885 """

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

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

888

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

890 """

891

892 # fd is closed by embracing read() routine

893 lines = fd.readlines()

894

895 atoms = None

896 cell = []

897 N = Nb = Nq = 0

898 scaled_positions = []

899 symbols = []

900 masses = []

901

902 # header

903 L = 0

904 while L < len(lines):

905

906 line = lines[L]

907

908 if 'Number of ions' in line:

909 N = int(line.split()[3])

910 elif 'Number of branches' in line:

911 Nb = int(line.split()[3])

912 elif 'Number of wavevectors' in line:

913 Nq = int(line.split()[3])

914 elif 'Unit cell vectors (A)' in line:

915 for _ in range(3):

916 L += 1

917 fields = lines[L].split()

918 cell.append([float(x) for x in fields[0:3]])

919 elif 'Fractional Co-ordinates' in line:

920 for _ in range(N):

921 L += 1

922 fields = lines[L].split()

923 scaled_positions.append([float(x) for x in fields[1:4]])

924 symbols.append(fields[4])

925 masses.append(float(fields[5]))

926 elif 'END header' in line:

927 L += 1

928 atoms = ase.Atoms(symbols=symbols,

929 scaled_positions=scaled_positions,

930 cell=cell)

931 break

932

933 L += 1

934

935 # Eigenmodes and -vectors

936 if frequency_factor is None:

937 Kayser_to_eV = 1E2 * 2 * np.pi * units['hbar'] * units['c']

938 # N.B. "fixed default" unit for frequencies in .phonon files is "cm-1"

939 # (i.e. the latter is unaffected by the internal unit conversion system of

940 # CASTEP!) set conversion factor to convert therefrom to eV by default for

941 # now

942 frequency_factor = Kayser_to_eV

943 qpoints = []

944 weights = []

945 frequencies = []

946 displacements = []

947 for _ in range(Nq):

948 fields = lines[L].split()

949 qpoints.append([float(x) for x in fields[2:5]])

950 weights.append(float(fields[5]))

951 freqs = []

952 for _ in range(Nb):

953 L += 1

954 fields = lines[L].split()

955 freqs.append(frequency_factor * float(fields[1]))

956 frequencies.append(np.array(freqs))

957

958 # skip the two Phonon Eigenvectors header lines

959 L += 2

960

961 # generate a list of displacements with a structure that is identical to

962 # what is stored internally in the Vibrations class (see in

963 # ase.vibrations.Vibrations.modes):

964 # np.array(displacements).shape == (Nb,3*N)

965

966 disps = []

967 for _ in range(Nb):

968 disp_coords = []

969 for _ in range(N):

970 L += 1

971 fields = lines[L].split()

972 disp_x = float(fields[2]) + float(fields[3]) * 1.0j

973 disp_y = float(fields[4]) + float(fields[5]) * 1.0j

974 disp_z = float(fields[6]) + float(fields[7]) * 1.0j

975 disp_coords.extend([disp_x, disp_y, disp_z])

976 disps.append(np.array(disp_coords))

977 displacements.append(np.array(disps))

978

979 if read_vib_data:

980 if gamma_only:

981 vibdata = [frequencies[0], displacements[0]]

982 else:

983 vibdata = [qpoints, weights, frequencies, displacements]

984 return vibdata, atoms

985 else:

986 return atoms

987

988

989def read_md(filename, index=None, return_scalars=False,

990 units=units_CODATA2002):

991 """Wrapper function for the more generic read() functionality.

992

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

994 only. For documentation see read_castep_md()

995 """

996 if return_scalars:

997 full_output = True

998 else:

999 full_output = False

1000

1001 from ase.io import read

1002 return read(filename, index=index, format='castep-md',

1003 full_output=full_output, return_scalars=return_scalars,

1004 units=units)

1005

1006

1007def read_castep_md(fd, index=None, return_scalars=False,

1008 units=units_CODATA2002):

1009 """Reads a .md file written by a CASTEP MolecularDynamics task

1010 and returns the trajectory stored therein as a list of atoms object.

1011

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

1013

1014 from ase.calculators.singlepoint import SinglePointCalculator

1015

1016 factors = {

1017 't': units['t0'] * 1E15, # fs

1018 'E': units['Eh'], # eV

1019 'T': units['Eh'] / units['kB'],

1020 'P': units['Eh'] / units['a0']**3 * units['Pascal'],

1021 'h': units['a0'],

1022 'hv': units['a0'] / units['t0'],

1023 'S': units['Eh'] / units['a0']**3,

1024 'R': units['a0'],

1025 'V': np.sqrt(units['Eh'] / units['me']),

1026 'F': units['Eh'] / units['a0']}

1027

1028 # fd is closed by embracing read() routine

1029 lines = fd.readlines()

1030

1031 L = 0

1032 while 'END header' not in lines[L]:

1033 L += 1

1034 l_end_header = L

1035 lines = lines[l_end_header + 1:]

1036 times = []

1037 energies = []

1038 temperatures = []

1039 pressures = []

1040 traj = []

1041

1042 # Initialization

1043 time = None

1044 Epot = None

1045 Ekin = None

1046 EH = None

1047 temperature = None

1048 pressure = None

1049 symbols = None

1050 positions = None

1051 cell = None

1052 velocities = None

1053 symbols = []

1054 positions = []

1055 velocities = []

1056 forces = []

1057 cell = np.eye(3)

1058 cell_velocities = []

1059 stress = []

1060

1061 for (L, line) in enumerate(lines):

1062 fields = line.split()

1063 if len(fields) == 0:

1064 if L != 0:

1065 times.append(time)

1066 energies.append([Epot, EH, Ekin])

1067 temperatures.append(temperature)

1068 pressures.append(pressure)

1069 atoms = ase.Atoms(symbols=symbols,

1070 positions=positions,

1071 cell=cell)

1072 atoms.set_velocities(velocities)

1073 if len(stress) == 0:

1074 atoms.calc = SinglePointCalculator(

1075 atoms=atoms, energy=Epot, forces=forces)

1076 else:

1077 atoms.calc = SinglePointCalculator(

1078 atoms=atoms, energy=Epot,

1079 forces=forces, stress=stress)

1080 traj.append(atoms)

1081 symbols = []

1082 positions = []

1083 velocities = []

1084 forces = []

1085 cell = []

1086 cell_velocities = []

1087 stress = []

1088 continue

1089 if len(fields) == 1:

1090 time = factors['t'] * float(fields[0])

1091 continue

1092

1093 if fields[-1] == 'E':

1094 E = [float(x) for x in fields[0:3]]

1095 Epot, EH, Ekin = (factors['E'] * Ei for Ei in E)

1096 continue

1097

1098 if fields[-1] == 'T':

1099 temperature = factors['T'] * float(fields[0])

1100 continue

1101

1102 # only printed in case of variable cell calculation or calculate_stress

1103 # explicitly requested

1104 if fields[-1] == 'P':

1105 pressure = factors['P'] * float(fields[0])

1106 continue

1107 if fields[-1] == 'h':

1108 h = [float(x) for x in fields[0:3]]

1109 cell.append([factors['h'] * hi for hi in h])

1110 continue

1111

1112 # only printed in case of variable cell calculation

1113 if fields[-1] == 'hv':

1114 hv = [float(x) for x in fields[0:3]]

1115 cell_velocities.append([factors['hv'] * hvi for hvi in hv])

1116 continue

1117

1118 # only printed in case of variable cell calculation

1119 if fields[-1] == 'S':

1120 S = [float(x) for x in fields[0:3]]

1121 stress.append([factors['S'] * Si for Si in S])

1122 continue

1123 if fields[-1] == 'R':

1124 symbols.append(fields[0])

1125 R = [float(x) for x in fields[2:5]]

1126 positions.append([factors['R'] * Ri for Ri in R])

1127 continue

1128 if fields[-1] == 'V':

1129 V = [float(x) for x in fields[2:5]]

1130 velocities.append([factors['V'] * Vi for Vi in V])

1131 continue

1132 if fields[-1] == 'F':

1133 F = [float(x) for x in fields[2:5]]

1134 forces.append([factors['F'] * Fi for Fi in F])

1135 continue

1136

1137 if index is None:

1138 pass

1139 else:

1140 traj = traj[index]

1141

1142 if return_scalars:

1143 data = [times, energies, temperatures, pressures]

1144 return data, traj

1145 else:

1146 return traj

1147

1148

1149# Routines that only the calculator requires

1150

1151def read_param(filename='', calc=None, fd=None, get_interface_options=False):

1152 if fd is None:

1153 if filename == '':

1154 raise ValueError('One between filename and fd must be provided')

1155 fd = open(filename)

1156 elif filename:

1157 warnings.warn('Filestream used to read param, file name will be '

1158 'ignored')

1159

1160 # If necessary, get the interface options

1161 if get_interface_options:

1162 int_opts = {}

1163 optre = re.compile(r'# ASE_INTERFACE ([^\s]+) : ([^\s]+)')

1164

1165 lines = fd.readlines()

1166 fd.seek(0)

1167

1168 for line in lines:

1169 m = optre.search(line)

1170 if m:

1171 int_opts[m.groups()[0]] = m.groups()[1]

1172

1173 data = read_freeform(fd)

1174

1175 if calc is None:

1176 from ase.calculators.castep import Castep

1177 calc = Castep(check_castep_version=False, keyword_tolerance=2)

1178

1179 for kw, (val, otype) in data.items():

1180 if otype == 'block':

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

1182 calc.param.__setattr__(kw, val)

1183

1184 if not get_interface_options:

1185 return calc

1186 else:

1187 return calc, int_opts

1188

1189

1190def write_param(filename, param, check_checkfile=False,

1191 force_write=False,

1192 interface_options=None):

1193 """Writes a CastepParam object to a CASTEP .param file

1194

1195 Parameters:

1196 filename: the location of the file to write to. If it

1197 exists it will be overwritten without warning. If it

1198 doesn't it will be created.

1199 param: a CastepParam instance

1200 check_checkfile : if set to True, write_param will

1201 only write continuation or reuse statement

1202 if a restart file exists in the same directory

1203 """

1204 if os.path.isfile(filename) and not force_write:

1205 warnings.warn('ase.io.castep.write_param: Set optional argument '

1206 'force_write=True to overwrite %s.' % filename)

1207 return False

1208

1209 out = paropen(filename, 'w')

1210 out.write('#######################################################\n')

1211 out.write(f'#CASTEP param file: {filename}\n')

1212 out.write('#Created using the Atomic Simulation Environment (ASE)#\n')

1213 if interface_options is not None:

1214 out.write('# Internal settings of the calculator\n')

1215 out.write('# This can be switched off by settings\n')

1216 out.write('# calc._export_settings = False\n')

1217 out.write('# If stated, this will be automatically processed\n')

1218 out.write('# by ase.io.castep.read_seed()\n')

1219 for option, value in sorted(interface_options.items()):

1220 out.write(f'# ASE_INTERFACE {option} : {value}\n')

1221 out.write('#######################################################\n\n')

1222

1223 if check_checkfile:

1224 param = deepcopy(param) # To avoid modifying the parent one

1225 for checktype in ['continuation', 'reuse']:

1226 opt = getattr(param, checktype)

1227 if opt and opt.value:

1228 fname = opt.value

1229 if fname == 'default':

1230 fname = os.path.splitext(filename)[0] + '.check'

1231 if not (os.path.exists(fname) or

1232 # CASTEP also understands relative path names, hence

1233 # also check relative to the param file directory

1234 os.path.exists(

1235 os.path.join(os.path.dirname(filename),

1236 opt.value))):

1237 opt.clear()

1238

1239 write_freeform(out, param)

1240

1241 out.close()

1242

1243

1244def read_seed(seed, new_seed=None, ignore_internal_keys=False):

1245 """A wrapper around the CASTEP Calculator in conjunction with

1246 read_cell and read_param. Basically this can be used to reuse

1247 a previous calculation which results in a triple of

1248 cell/param/castep file. The label of the calculation if pre-

1249 fixed with `copy_of_` and everything else will be recycled as

1250 much as possible from the addressed calculation.

1251

1252 Please note that this routine will return an atoms ordering as specified

1253 in the cell file! It will thus undo the potential reordering internally

1254 done by castep.

1255 """

1256

1257 directory = os.path.abspath(os.path.dirname(seed))

1258 seed = os.path.basename(seed)

1259

1260 paramfile = os.path.join(directory, f'{seed}.param')

1261 cellfile = os.path.join(directory, f'{seed}.cell')

1262 castepfile = os.path.join(directory, f'{seed}.castep')

1263 checkfile = os.path.join(directory, f'{seed}.check')

1264

1265 atoms = read_cell(cellfile)

1266 atoms.calc._directory = directory

1267 atoms.calc._rename_existing_dir = False

1268 atoms.calc._castep_pp_path = directory

1269 atoms.calc.merge_param(paramfile,

1270 ignore_internal_keys=ignore_internal_keys)

1271 if new_seed is None:

1272 atoms.calc._label = f'copy_of_{seed}'

1273 else:

1274 atoms.calc._label = str(new_seed)

1275 if os.path.isfile(castepfile):

1276 # _set_atoms needs to be True here

1277 # but we set it right back to False

1278 # atoms.calc._set_atoms = False

1279 # BUGFIX: I do not see a reason to do that!

1280 atoms.calc.read(castepfile)

1281 # atoms.calc._set_atoms = False

1282

1283 # if here is a check file, we also want to re-use this information

1284 if os.path.isfile(checkfile):

1285 atoms.calc._check_file = os.path.basename(checkfile)

1286

1287 # sync the top-level object with the

1288 # one attached to the calculator

1289 atoms = atoms.calc.atoms

1290 else:

1291 # There are cases where we only want to restore a calculator/atoms

1292 # setting without a castep file...

1293 # No print statement required in these cases

1294 warnings.warn(

1295 'Corresponding *.castep file not found. '

1296 'Atoms object will be restored from *.cell and *.param only.')

1297 atoms.calc.push_oldstate()

1298

1299 return atoms

1300

1301

1302@reader

1303def read_bands(fd, units=units_CODATA2002):

1304 """Read Castep.bands file to kpoints, weights and eigenvalues.

1305

1306 Parameters

1307 ----------

1308 fd : str | io.TextIOBase

1309 Path to the `.bands` file or file descriptor for open `.bands` file.

1310 units : dict

1311 Conversion factors for atomic units.

1312

1313 Returns

1314 -------

1315 kpts : np.ndarray

1316 1d NumPy array for k-point coordinates.

1317 weights : np.ndarray

1318 1d NumPy array for k-point weights.

1319 eigenvalues : np.ndarray

1320 NumPy array for eigenvalues with shape (spin, kpts, nbands).

1321 efermi : float

1322 Fermi energy.

1323

1324 """

1325 Hartree = units['Eh']

1326

1327 nkpts = int(fd.readline().split()[-1])

1328 nspin = int(fd.readline().split()[-1])

1329 _ = float(fd.readline().split()[-1])

1330 nbands = int(fd.readline().split()[-1])

1331 efermi = float(fd.readline().split()[-1])

1332

1333 kpts = np.zeros((nkpts, 3))

1334 weights = np.zeros(nkpts)

1335 eigenvalues = np.zeros((nspin, nkpts, nbands))

1336

1337 # Skip unit cell

1338 for _ in range(4):

1339 fd.readline()

1340

1341 def _kptline_to_i_k_wt(line: str) -> Tuple[int, List[float], float]:

1342 split_line = line.split()

1343 i_kpt = int(split_line[1]) - 1

1344 kpt = list(map(float, split_line[2:5]))

1345 wt = float(split_line[5])

1346 return i_kpt, kpt, wt

1347

1348 # CASTEP often writes these out-of-order, so check index and write directly

1349 # to the correct row

1350 for _ in range(nkpts):

1351 i_kpt, kpt, wt = _kptline_to_i_k_wt(fd.readline())

1352 kpts[i_kpt, :], weights[i_kpt] = kpt, wt

1353 for spin in range(nspin):

1354 fd.readline() # Skip 'Spin component N' line

1355 eigenvalues[spin, i_kpt, :] = [float(fd.readline())

1356 for _ in range(nbands)]

1357

1358 return (kpts, weights, eigenvalues * Hartree, efermi * Hartree)

Coverage for /builds/alexhroom/ase/ase/io/castep/init.py: 50.22%

695 statements