Coverage for /builds/hweiske/ase/ase/io/castep/__init__.py: 43.35%
752 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +0000
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.parallel import paropen
25from ase.spacegroup import Spacegroup
26from ase.utils import atoms_to_spglib_cell
28units_ase = {
29 'hbar': ase.units._hbar * ase.units.J,
30 'Eh': ase.units.Hartree,
31 'kB': ase.units.kB,
32 'a0': ase.units.Bohr,
33 't0': ase.units._hbar * ase.units.J / ase.units.Hartree,
34 'c': ase.units._c,
35 'me': ase.units._me / ase.units._amu,
36 'Pascal': 1.0 / ase.units.Pascal}
38# CODATA1986 (included herein for the sake of completeness)
39# taken from
40# http://physics.nist.gov/cuu/Archive/1986RMP.pdf
41units_CODATA1986 = {
42 'hbar': 6.5821220E-16, # eVs
43 'Eh': 27.2113961, # eV
44 'kB': 8.617385E-5, # eV/K
45 'a0': 0.529177249, # A
46 'c': 299792458, # m/s
47 'e': 1.60217733E-19, # C
48 'me': 5.485799110E-4} # u
50# CODATA2002: default in CASTEP 5.01
51# (-> check in more recent CASTEP in case of numerical discrepancies?!)
52# taken from
53# http://physics.nist.gov/cuu/Document/all_2002.pdf
54units_CODATA2002 = {
55 'hbar': 6.58211915E-16, # eVs
56 'Eh': 27.2113845, # eV
57 'kB': 8.617343E-5, # eV/K
58 'a0': 0.5291772108, # A
59 'c': 299792458, # m/s
60 'e': 1.60217653E-19, # C
61 'me': 5.4857990945E-4} # u
63# (common) derived entries
64for d in (units_CODATA1986, units_CODATA2002):
65 d['t0'] = d['hbar'] / d['Eh'] # s
66 d['Pascal'] = d['e'] * 1E30 # Pa
69__all__ = [
70 # routines for the generic io function
71 'read_castep',
72 'read_castep_castep',
73 'read_castep_castep_old',
74 'read_cell',
75 'read_castep_cell',
76 'read_geom',
77 'read_castep_geom',
78 'read_phonon',
79 'read_castep_phonon',
80 # additional reads that still need to be wrapped
81 'read_md',
82 'read_param',
83 'read_seed',
84 # write that is already wrapped
85 'write_castep_cell',
86 # param write - in principle only necessary in junction with the calculator
87 'write_param']
90def write_freeform(fd, outputobj):
91 """
92 Prints out to a given file a CastepInputFile or derived class, such as
93 CastepCell or CastepParam.
94 """
96 options = outputobj._options
98 # Some keywords, if present, are printed in this order
99 preferred_order = ['lattice_cart', 'lattice_abc',
100 'positions_frac', 'positions_abs',
101 'species_pot', 'symmetry_ops', # CELL file
102 'task', 'cut_off_energy' # PARAM file
103 ]
105 keys = outputobj.get_attr_dict().keys()
106 # This sorts only the ones in preferred_order and leaves the rest
107 # untouched
108 keys = sorted(keys, key=lambda x: preferred_order.index(x)
109 if x in preferred_order
110 else len(preferred_order))
112 for kw in keys:
113 opt = options[kw]
114 if opt.type.lower() == 'block':
115 fd.write('%BLOCK {0}\n{1}\n%ENDBLOCK {0}\n\n'.format(
116 kw.upper(),
117 opt.value.strip('\n')))
118 else:
119 fd.write(f'{kw.upper()}: {opt.value}\n')
122def write_cell(filename, atoms, positions_frac=False, castep_cell=None,
123 force_write=False):
124 """
125 Wrapper function for the more generic write() functionality.
127 Note that this is function is intended to maintain backwards-compatibility
128 only.
129 """
130 from ase.io import write
132 write(filename, atoms, positions_frac=positions_frac,
133 castep_cell=castep_cell, force_write=force_write)
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.
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
148 Note that the "force_write" keywords has no effect currently.
150 Arguments:
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 """
166 if isinstance(atoms, list):
167 if len(atoms) > 1:
168 atoms = atoms[-1]
170 # Header
171 fd.write('# written by ASE\n\n')
173 # To write this we simply use the existing Castep calculator, or create
174 # one
175 from ase.calculators.castep import Castep, CastepCell
177 try:
178 has_cell = isinstance(atoms.calc.cell, CastepCell)
179 except AttributeError:
180 has_cell = False
182 if has_cell:
183 cell = deepcopy(atoms.calc.cell)
184 else:
185 cell = Castep(keyword_tolerance=2).cell
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()]
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()
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'):
206 from ase.data import atomic_masses
207 masses = atoms.get_array('masses')
208 custom_masses = {}
210 for i, species in enumerate(elems):
211 custom_mass = masses[i]
213 # build record of different masses for each species
214 if species not in custom_masses:
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]}
221 # if multiple masses found for a species
222 elif custom_mass not in custom_masses[species].keys():
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]))
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)
247 # create species_mass block
248 mass_block = []
250 for el, mass_dict in custom_masses.items():
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))
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))
279 cell.species_mass = mass_block
281 if atoms.has('castep_labels'):
282 labels = atoms.get_array('castep_labels')
283 else:
284 labels = ['NULL'] * len(elems)
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)
295 pos_block = []
296 pos_block_format = '%s ' + cell_block_format
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)
308 setattr(cell, pos_keyword, pos_block)
310 constr_block = _make_block_ionic_constraints(atoms)
311 if constr_block:
312 cell.ionic_constraints = constr_block
314 write_freeform(fd, cell)
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
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
365def _calc_normal_vectors(constr: FixedLine) -> Tuple[np.ndarray, np.ndarray]:
366 direction = constr.dir
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)
376 n2 = np.cross(direction, n1)
377 n2 = n2 / np.linalg.norm(n2)
379 return n1, n2
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 """
389 from ase.io.castep.castep_input_file import CastepInputFile
391 inputobj = CastepInputFile(keyword_tolerance=2)
393 filelines = fd.readlines()
395 keyw = None
396 read_block = False
397 block_lines = None
399 for i, l in enumerate(filelines):
401 # Strip all comments, aka anything after a hash
402 L = re.split(r'[#!;]', l, 1)[0].strip()
404 if L == '':
405 # Empty line... skip
406 continue
408 lsplit = re.split(r'\s*[:=]*\s+', L, 1)
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
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()
434 # Now save the value
435 if read_block:
436 block_lines = []
437 else:
438 inputobj.__setattr__(keyw, ' '.join(lsplit[1:]))
440 return inputobj.get_attr_dict(types=True)
443def read_cell(filename, index=None):
444 """
445 Wrapper function for the more generic read() functionality.
447 Note that this is function is intended to maintain backwards-compatibility
448 only.
449 """
450 from ase.io import read
451 return read(filename, index=index, format='castep-cell')
454def read_castep_cell(fd, index=None, calculator_args={}, find_spg=False,
455 units=units_CODATA2002):
456 """Read a .cell file and return an atoms object.
457 Any value found that does not fit the atoms API
458 will be stored in the atoms.calc attribute.
460 By default, the Castep calculator will be tolerant and in the absence of a
461 castep_keywords.json file it will just accept all keywords that aren't
462 automatically parsed.
463 """
465 from ase.calculators.castep import Castep
467 cell_units = { # Units specifiers for CASTEP
468 'bohr': units_CODATA2002['a0'],
469 'ang': 1.0,
470 'm': 1e10,
471 'cm': 1e8,
472 'nm': 10,
473 'pm': 1e-2
474 }
476 calc = Castep(**calculator_args)
478 if calc.cell.castep_version == 0 and calc._kw_tol < 3:
479 # No valid castep_keywords.json was found
480 warnings.warn(
481 'read_cell: Warning - Was not able to validate CASTEP input. '
482 'This may be due to a non-existing '
483 '"castep_keywords.json" '
484 'file or a non-existing CASTEP installation. '
485 'Parsing will go on but keywords will not be '
486 'validated and may cause problems if incorrect during a CASTEP '
487 'run.')
489 celldict = read_freeform(fd)
491 def parse_blockunit(line_tokens, blockname):
492 u = 1.0
493 if len(line_tokens[0]) == 1:
494 usymb = line_tokens[0][0].lower()
495 u = cell_units.get(usymb, 1)
496 if usymb not in cell_units:
497 warnings.warn('read_cell: Warning - ignoring invalid '
498 'unit specifier in %BLOCK {} '
499 '(assuming Angstrom instead)'.format(blockname))
500 line_tokens = line_tokens[1:]
501 return u, line_tokens
503 # Arguments to pass to the Atoms object at the end
504 aargs = {
505 'pbc': True
506 }
508 # Start by looking for the lattice
509 lat_keywords = [w in celldict for w in ('lattice_cart', 'lattice_abc')]
510 if all(lat_keywords):
511 warnings.warn('read_cell: Warning - two lattice blocks present in the'
512 ' same file. LATTICE_ABC will be ignored')
513 elif not any(lat_keywords):
514 raise ValueError('Cell file must contain at least one between '
515 'LATTICE_ABC and LATTICE_CART')
517 if 'lattice_abc' in celldict:
519 lines = celldict.pop('lattice_abc')[0].split('\n')
520 line_tokens = [line.split() for line in lines]
522 u, line_tokens = parse_blockunit(line_tokens, 'lattice_abc')
524 if len(line_tokens) != 2:
525 warnings.warn('read_cell: Warning - ignoring additional '
526 'lines in invalid %BLOCK LATTICE_ABC')
528 abc = [float(p) * u for p in line_tokens[0][:3]]
529 angles = [float(phi) for phi in line_tokens[1][:3]]
531 aargs['cell'] = cellpar_to_cell(abc + angles)
533 if 'lattice_cart' in celldict:
535 lines = celldict.pop('lattice_cart')[0].split('\n')
536 line_tokens = [line.split() for line in lines]
538 u, line_tokens = parse_blockunit(line_tokens, 'lattice_cart')
540 if len(line_tokens) != 3:
541 warnings.warn('read_cell: Warning - ignoring more than '
542 'three lattice vectors in invalid %BLOCK '
543 'LATTICE_CART')
545 aargs['cell'] = [[float(x) * u for x in lt[:3]] for lt in line_tokens]
547 # Now move on to the positions
548 pos_keywords = [w in celldict
549 for w in ('positions_abs', 'positions_frac')]
551 if all(pos_keywords):
552 warnings.warn('read_cell: Warning - two lattice blocks present in the'
553 ' same file. POSITIONS_FRAC will be ignored')
554 del celldict['positions_frac']
555 elif not any(pos_keywords):
556 raise ValueError('Cell file must contain at least one between '
557 'POSITIONS_FRAC and POSITIONS_ABS')
559 aargs['symbols'] = []
560 pos_type = 'positions'
561 pos_block = celldict.pop('positions_abs', [None])[0]
562 if pos_block is None:
563 pos_type = 'scaled_positions'
564 pos_block = celldict.pop('positions_frac', [None])[0]
565 aargs[pos_type] = []
567 lines = pos_block.split('\n')
568 line_tokens = [line.split() for line in lines]
570 if 'scaled' not in pos_type:
571 u, line_tokens = parse_blockunit(line_tokens, 'positions_abs')
572 else:
573 u = 1.0
575 # Here we extract all the possible additional info
576 # These are marked by their type
578 add_info = {
579 'SPIN': (float, 0.0), # (type, default)
580 'MAGMOM': (float, 0.0),
581 'LABEL': (str, 'NULL')
582 }
583 add_info_arrays = {k: [] for k in add_info}
585 def parse_info(raw_info):
587 re_keys = (r'({0})\s*[=:\s]{{1}}\s'
588 r'*([^\s]*)').format('|'.join(add_info.keys()))
589 # Capture all info groups
590 info = re.findall(re_keys, raw_info)
591 info = {g[0]: add_info[g[0]][0](g[1]) for g in info}
592 return info
594 # Array for custom species (a CASTEP special thing)
595 # Usually left unused
596 custom_species = None
598 for tokens in line_tokens:
599 # Now, process the whole 'species' thing
600 spec_custom = tokens[0].split(':', 1)
601 elem = spec_custom[0]
602 if len(spec_custom) > 1 and custom_species is None:
603 # Add it to the custom info!
604 custom_species = list(aargs['symbols'])
605 if custom_species is not None:
606 custom_species.append(tokens[0])
607 aargs['symbols'].append(elem)
608 aargs[pos_type].append([float(p) * u for p in tokens[1:4]])
609 # Now for the additional information
610 info = ' '.join(tokens[4:])
611 info = parse_info(info)
612 for k in add_info:
613 add_info_arrays[k] += [info.get(k, add_info[k][1])]
615 # read in custom species mass
616 if 'species_mass' in celldict:
617 spec_list = custom_species if custom_species else aargs['symbols']
618 aargs['masses'] = [None for _ in spec_list]
619 lines = celldict.pop('species_mass')[0].split('\n')
620 line_tokens = [line.split() for line in lines]
622 if len(line_tokens[0]) == 1:
623 if line_tokens[0][0].lower() not in ('amu', 'u'):
624 raise ValueError(
625 "unit specifier '{}' in %BLOCK SPECIES_MASS "
626 "not recognised".format(
627 line_tokens[0][0].lower()))
628 line_tokens = line_tokens[1:]
630 for tokens in line_tokens:
631 token_pos_list = [i for i, x in enumerate(
632 spec_list) if x == tokens[0]]
633 if len(token_pos_list) == 0:
634 warnings.warn(
635 'read_cell: Warning - ignoring unused '
636 'species mass {} in %BLOCK SPECIES_MASS'.format(
637 tokens[0]))
638 for idx in token_pos_list:
639 aargs['masses'][idx] = tokens[1]
641 # Now on to the species potentials...
642 if 'species_pot' in celldict:
643 lines = celldict.pop('species_pot')[0].split('\n')
644 line_tokens = [line.split() for line in lines]
646 for tokens in line_tokens:
647 if len(tokens) == 1:
648 # It's a library
649 all_spec = (set(custom_species) if custom_species is not None
650 else set(aargs['symbols']))
651 for s in all_spec:
652 calc.cell.species_pot = (s, tokens[0])
653 else:
654 calc.cell.species_pot = tuple(tokens[:2])
656 # Ionic constraints
657 raw_constraints = {}
659 if 'ionic_constraints' in celldict:
660 lines = celldict.pop('ionic_constraints')[0].split('\n')
661 line_tokens = [line.split() for line in lines]
663 for tokens in line_tokens:
664 if len(tokens) != 6:
665 continue
666 _, species, nic, x, y, z = tokens
667 # convert xyz to floats
668 x = float(x)
669 y = float(y)
670 z = float(z)
672 nic = int(nic)
673 if (species, nic) not in raw_constraints:
674 raw_constraints[(species, nic)] = []
675 raw_constraints[(species, nic)].append(np.array(
676 [x, y, z]))
678 # Symmetry operations
679 if 'symmetry_ops' in celldict:
680 lines = celldict.pop('symmetry_ops')[0].split('\n')
681 line_tokens = [line.split() for line in lines]
683 # Read them in blocks of four
684 blocks = np.array(line_tokens).astype(float)
685 if (len(blocks.shape) != 2 or blocks.shape[1] != 3
686 or blocks.shape[0] % 4 != 0):
687 warnings.warn('Warning: could not parse SYMMETRY_OPS'
688 ' block properly, skipping')
689 else:
690 blocks = blocks.reshape((-1, 4, 3))
691 rotations = blocks[:, :3]
692 translations = blocks[:, 3]
694 # Regardless of whether we recognize them, store these
695 calc.cell.symmetry_ops = (rotations, translations)
697 # Anything else that remains, just add it to the cell object:
698 for k, (val, otype) in celldict.items():
699 try:
700 if otype == 'block':
701 val = val.split('\n') # Avoids a bug for one-line blocks
702 calc.cell.__setattr__(k, val)
703 except Exception as e:
704 raise RuntimeError(
705 f'Problem setting calc.cell.{k} = {val}: {e}')
707 # Get the relevant additional info
708 aargs['magmoms'] = np.array(add_info_arrays['SPIN'])
709 # SPIN or MAGMOM are alternative keywords
710 aargs['magmoms'] = np.where(aargs['magmoms'] != 0,
711 aargs['magmoms'],
712 add_info_arrays['MAGMOM'])
713 labels = np.array(add_info_arrays['LABEL'])
715 aargs['calculator'] = calc
717 atoms = ase.Atoms(**aargs)
719 # Spacegroup...
720 if find_spg:
721 # Try importing spglib
722 try:
723 import spglib
724 except ImportError:
725 warnings.warn('spglib not found installed on this system - '
726 'automatic spacegroup detection is not possible')
727 spglib = None
729 if spglib is not None:
730 symmd = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms))
731 atoms_spg = Spacegroup(int(symmd['number']))
732 atoms.info['spacegroup'] = atoms_spg
734 atoms.new_array('castep_labels', labels)
735 if custom_species is not None:
736 atoms.new_array('castep_custom_species', np.array(custom_species))
738 fixed_atoms = []
739 constraints = []
740 index_dict = atoms.symbols.indices()
741 for (species, nic), value in raw_constraints.items():
743 absolute_nr = index_dict[species][nic - 1]
744 if len(value) == 3:
745 # Check if they are linearly independent
746 if np.linalg.det(value) == 0:
747 warnings.warn(
748 'Error: Found linearly dependent constraints attached '
749 'to atoms %s' %
750 (absolute_nr))
751 continue
752 fixed_atoms.append(absolute_nr)
753 elif len(value) == 2:
754 direction = np.cross(value[0], value[1])
755 # Check if they are linearly independent
756 if np.linalg.norm(direction) == 0:
757 warnings.warn(
758 'Error: Found linearly dependent constraints attached '
759 'to atoms %s' %
760 (absolute_nr))
761 continue
762 constraint = FixedLine(indices=absolute_nr, direction=direction)
763 constraints.append(constraint)
764 elif len(value) == 1:
765 direction = np.array(value[0], dtype=float)
766 constraint = FixedPlane(indices=absolute_nr, direction=direction)
767 constraints.append(constraint)
768 else:
769 warnings.warn(
770 f'Error: Found {len(value)} statements attached to atoms '
771 f'{absolute_nr}'
772 )
774 # we need to sort the fixed atoms list in order not to raise an assertion
775 # error in FixAtoms
776 if fixed_atoms:
777 constraints.append(FixAtoms(indices=sorted(fixed_atoms)))
778 if constraints:
779 atoms.set_constraint(constraints)
781 atoms.calc.atoms = atoms
782 atoms.calc.push_oldstate()
784 return atoms
787def read_castep(filename, index=None):
788 """
789 Wrapper function for the more generic read() functionality.
791 Note that this is function is intended to maintain backwards-compatibility
792 only.
793 """
794 from ase.io import read
795 return read(filename, index=index, format='castep-castep')
798def read_castep_castep(fd, index=None):
799 """
800 Reads a .castep file and returns an atoms object.
801 The calculator information will be stored in the calc attribute.
803 There is no use of the "index" argument as of now, it is just inserted for
804 convenience to comply with the generic "read()" in ase.io
806 Please note that this routine will return an atom ordering as found
807 within the castep file. This means that the species will be ordered by
808 ascending atomic numbers. The atoms witin a species are ordered as given
809 in the original cell file.
811 Note: This routine returns a single atoms_object only, the last
812 configuration in the file. Yet, if you want to parse an MD run, use the
813 novel function `read_md()`
814 """
816 from ase.calculators.castep import Castep
818 try:
819 calc = Castep()
820 except Exception as e:
821 # No CASTEP keywords found?
822 warnings.warn(f'WARNING: {e} Using fallback .castep reader...')
823 # Fall back on the old method
824 return read_castep_castep_old(fd, index)
826 calc.read(castep_file=fd)
828 # now we trick the calculator instance such that we can savely extract
829 # energies and forces from this atom. Basically what we do is to trick the
830 # internal routine calculation_required() to always return False such that
831 # we do not need to re-run a CASTEP calculation.
832 #
833 # Probably we can solve this with a flag to the read() routine at some
834 # point, but for the moment I do not want to change too much in there.
835 calc._old_atoms = calc.atoms
836 calc._old_param = calc.param
837 calc._old_cell = calc.cell
839 return [calc.atoms] # Returning in the form of a list for next()
842def read_castep_castep_old(fd, index=None):
843 """
844 DEPRECATED
845 Now replaced by ase.calculators.castep.Castep.read(). Left in for future
846 reference and backwards compatibility needs, as well as a fallback for
847 when castep_keywords.py can't be created.
849 Reads a .castep file and returns an atoms object.
850 The calculator information will be stored in the calc attribute.
851 If more than one SCF step is found, a list of all steps
852 will be stored in the traj attribute.
854 Note that the index argument has no effect as of now.
856 Please note that this routine will return an atom ordering as found
857 within the castep file. This means that the species will be ordered by
858 ascending atomic numbers. The atoms witin a species are ordered as given
859 in the original cell file.
860 """
861 from ase.calculators.singlepoint import SinglePointCalculator
863 lines = fd.readlines()
865 traj = []
866 energy_total = None
867 energy_0K = None
868 for i, line in enumerate(lines):
869 if 'NB est. 0K energy' in line:
870 energy_0K = float(line.split()[6])
871 # support also for dispersion correction
872 elif 'NB dispersion corrected est. 0K energy*' in line:
873 energy_0K = float(line.split()[-2])
874 elif 'Final energy, E' in line:
875 energy_total = float(line.split()[4])
876 elif 'Dispersion corrected final energy' in line:
877 pass
878 # dispcorr_energy_total = float(line.split()[-2])
879 # sedc_apply = True
880 elif 'Dispersion corrected final free energy' in line:
881 pass # dispcorr_energy_free = float(line.split()[-2])
882 elif 'dispersion corrected est. 0K energy' in line:
883 pass # dispcorr_energy_0K = float(line.split()[-2])
884 elif 'Unit Cell' in line:
885 cell = [x.split()[0:3] for x in lines[i + 3:i + 6]]
886 cell = np.array([[float(col) for col in row] for row in cell])
887 elif 'Cell Contents' in line:
888 geom_starts = i
889 start_found = False
890 for j, jline in enumerate(lines[geom_starts:]):
891 if jline.find('xxxxx') > 0 and start_found:
892 geom_stop = j + geom_starts
893 break
894 if jline.find('xxxx') > 0 and not start_found:
895 geom_start = j + geom_starts + 4
896 start_found = True
897 species = [line.split()[1] for line in lines[geom_start:geom_stop]]
898 geom = np.dot(np.array([[float(col) for col in line.split()[3:6]]
899 for line in lines[geom_start:geom_stop]]),
900 cell)
901 elif 'Writing model to' in line:
902 atoms = ase.Atoms(
903 cell=cell,
904 pbc=True,
905 positions=geom,
906 symbols=''.join(species))
907 # take 0K energy where available, else total energy
908 if energy_0K:
909 energy = energy_0K
910 else:
911 energy = energy_total
912 # generate a minimal single-point calculator
913 sp_calc = SinglePointCalculator(atoms=atoms,
914 energy=energy,
915 forces=None,
916 magmoms=None,
917 stress=None)
918 atoms.calc = sp_calc
919 traj.append(atoms)
920 if index is None:
921 return traj
922 else:
923 return traj[index]
926def read_geom(filename, index=':', units=units_CODATA2002):
927 """
928 Wrapper function for the more generic read() functionality.
930 Note that this is function is intended to maintain backwards-compatibility
931 only. Keyword arguments will be passed to read_castep_geom().
932 """
933 from ase.io import read
934 return read(filename, index=index, format='castep-geom', units=units)
937def read_castep_geom(fd, index=None, units=units_CODATA2002):
938 """Reads a .geom file produced by the CASTEP GeometryOptimization task and
939 returns an atoms object.
940 The information about total free energy and forces of each atom for every
941 relaxation step will be stored for further analysis especially in a
942 single-point calculator.
943 Note that everything in the .geom file is in atomic units, which has
944 been conversed to commonly used unit angstrom(length) and eV (energy).
946 Note that the index argument has no effect as of now.
948 Contribution by Wei-Bing Zhang. Thanks!
950 Routine now accepts a filedescriptor in order to out-source the gz and
951 bz2 handling to formats.py. Note that there is a fallback routine
952 read_geom() that behaves like previous versions did.
953 """
954 from ase.calculators.singlepoint import SinglePointCalculator
956 # fd is closed by embracing read() routine
957 txt = fd.readlines()
959 traj = []
961 Hartree = units['Eh']
962 Bohr = units['a0']
964 # Yeah, we know that...
965 # print('N.B.: Energy in .geom file is not 0K extrapolated.')
966 for i, line in enumerate(txt):
967 if line.find('<-- E') > 0:
968 start_found = True
969 energy = float(line.split()[0]) * Hartree
970 cell = [x.split()[0:3] for x in txt[i + 1:i + 4]]
971 cell = np.array([[float(col) * Bohr for col in row] for row in
972 cell])
973 if line.find('<-- R') > 0 and start_found:
974 start_found = False
975 geom_start = i
976 for i, line in enumerate(txt[geom_start:]):
977 if line.find('<-- F') > 0:
978 geom_stop = i + geom_start
979 break
980 species = [line.split()[0] for line in
981 txt[geom_start:geom_stop]]
982 geom = np.array([[float(col) * Bohr for col in
983 line.split()[2:5]] for line in
984 txt[geom_start:geom_stop]])
985 forces = np.array([[float(col) * Hartree / Bohr for col in
986 line.split()[2:5]] for line in
987 txt[geom_stop:geom_stop
988 + (geom_stop - geom_start)]])
989 image = ase.Atoms(species, geom, cell=cell, pbc=True)
990 image.calc = SinglePointCalculator(
991 atoms=image, energy=energy, forces=forces)
992 traj.append(image)
994 if index is None:
995 return traj
996 else:
997 return traj[index]
1000def read_phonon(filename, index=None, read_vib_data=False,
1001 gamma_only=True, frequency_factor=None,
1002 units=units_CODATA2002):
1003 """
1004 Wrapper function for the more generic read() functionality.
1006 Note that this is function is intended to maintain backwards-compatibility
1007 only. For documentation see read_castep_phonon().
1008 """
1009 from ase.io import read
1011 if read_vib_data:
1012 full_output = True
1013 else:
1014 full_output = False
1016 return read(filename, index=index, format='castep-phonon',
1017 full_output=full_output, read_vib_data=read_vib_data,
1018 gamma_only=gamma_only, frequency_factor=frequency_factor,
1019 units=units)
1022def read_castep_phonon(fd, index=None, read_vib_data=False,
1023 gamma_only=True, frequency_factor=None,
1024 units=units_CODATA2002):
1025 """
1026 Reads a .phonon file written by a CASTEP Phonon task and returns an atoms
1027 object, as well as the calculated vibrational data if requested.
1029 Note that the index argument has no effect as of now.
1030 """
1032 # fd is closed by embracing read() routine
1033 lines = fd.readlines()
1035 atoms = None
1036 cell = []
1037 N = Nb = Nq = 0
1038 scaled_positions = []
1039 symbols = []
1040 masses = []
1042 # header
1043 L = 0
1044 while L < len(lines):
1046 line = lines[L]
1048 if 'Number of ions' in line:
1049 N = int(line.split()[3])
1050 elif 'Number of branches' in line:
1051 Nb = int(line.split()[3])
1052 elif 'Number of wavevectors' in line:
1053 Nq = int(line.split()[3])
1054 elif 'Unit cell vectors (A)' in line:
1055 for _ in range(3):
1056 L += 1
1057 fields = lines[L].split()
1058 cell.append([float(x) for x in fields[0:3]])
1059 elif 'Fractional Co-ordinates' in line:
1060 for _ in range(N):
1061 L += 1
1062 fields = lines[L].split()
1063 scaled_positions.append([float(x) for x in fields[1:4]])
1064 symbols.append(fields[4])
1065 masses.append(float(fields[5]))
1066 elif 'END header' in line:
1067 L += 1
1068 atoms = ase.Atoms(symbols=symbols,
1069 scaled_positions=scaled_positions,
1070 cell=cell)
1071 break
1073 L += 1
1075 # Eigenmodes and -vectors
1076 if frequency_factor is None:
1077 Kayser_to_eV = 1E2 * 2 * np.pi * units['hbar'] * units['c']
1078 # N.B. "fixed default" unit for frequencies in .phonon files is "cm-1"
1079 # (i.e. the latter is unaffected by the internal unit conversion system of
1080 # CASTEP!) set conversion factor to convert therefrom to eV by default for
1081 # now
1082 frequency_factor = Kayser_to_eV
1083 qpoints = []
1084 weights = []
1085 frequencies = []
1086 displacements = []
1087 for _ in range(Nq):
1088 fields = lines[L].split()
1089 qpoints.append([float(x) for x in fields[2:5]])
1090 weights.append(float(fields[5]))
1091 freqs = []
1092 for _ in range(Nb):
1093 L += 1
1094 fields = lines[L].split()
1095 freqs.append(frequency_factor * float(fields[1]))
1096 frequencies.append(np.array(freqs))
1098 # skip the two Phonon Eigenvectors header lines
1099 L += 2
1101 # generate a list of displacements with a structure that is identical to
1102 # what is stored internally in the Vibrations class (see in
1103 # ase.vibrations.Vibrations.modes):
1104 # np.array(displacements).shape == (Nb,3*N)
1106 disps = []
1107 for _ in range(Nb):
1108 disp_coords = []
1109 for _ in range(N):
1110 L += 1
1111 fields = lines[L].split()
1112 disp_x = float(fields[2]) + float(fields[3]) * 1.0j
1113 disp_y = float(fields[4]) + float(fields[5]) * 1.0j
1114 disp_z = float(fields[6]) + float(fields[7]) * 1.0j
1115 disp_coords.extend([disp_x, disp_y, disp_z])
1116 disps.append(np.array(disp_coords))
1117 displacements.append(np.array(disps))
1119 if read_vib_data:
1120 if gamma_only:
1121 vibdata = [frequencies[0], displacements[0]]
1122 else:
1123 vibdata = [qpoints, weights, frequencies, displacements]
1124 return vibdata, atoms
1125 else:
1126 return atoms
1129def read_md(filename, index=None, return_scalars=False,
1130 units=units_CODATA2002):
1131 """Wrapper function for the more generic read() functionality.
1133 Note that this function is intended to maintain backwards-compatibility
1134 only. For documentation see read_castep_md()
1135 """
1136 if return_scalars:
1137 full_output = True
1138 else:
1139 full_output = False
1141 from ase.io import read
1142 return read(filename, index=index, format='castep-md',
1143 full_output=full_output, return_scalars=return_scalars,
1144 units=units)
1147def read_castep_md(fd, index=None, return_scalars=False,
1148 units=units_CODATA2002):
1149 """Reads a .md file written by a CASTEP MolecularDynamics task
1150 and returns the trajectory stored therein as a list of atoms object.
1152 Note that the index argument has no effect as of now."""
1154 from ase.calculators.singlepoint import SinglePointCalculator
1156 factors = {
1157 't': units['t0'] * 1E15, # fs
1158 'E': units['Eh'], # eV
1159 'T': units['Eh'] / units['kB'],
1160 'P': units['Eh'] / units['a0']**3 * units['Pascal'],
1161 'h': units['a0'],
1162 'hv': units['a0'] / units['t0'],
1163 'S': units['Eh'] / units['a0']**3,
1164 'R': units['a0'],
1165 'V': np.sqrt(units['Eh'] / units['me']),
1166 'F': units['Eh'] / units['a0']}
1168 # fd is closed by embracing read() routine
1169 lines = fd.readlines()
1171 L = 0
1172 while 'END header' not in lines[L]:
1173 L += 1
1174 l_end_header = L
1175 lines = lines[l_end_header + 1:]
1176 times = []
1177 energies = []
1178 temperatures = []
1179 pressures = []
1180 traj = []
1182 # Initialization
1183 time = None
1184 Epot = None
1185 Ekin = None
1186 EH = None
1187 temperature = None
1188 pressure = None
1189 symbols = None
1190 positions = None
1191 cell = None
1192 velocities = None
1193 symbols = []
1194 positions = []
1195 velocities = []
1196 forces = []
1197 cell = np.eye(3)
1198 cell_velocities = []
1199 stress = []
1201 for (L, line) in enumerate(lines):
1202 fields = line.split()
1203 if len(fields) == 0:
1204 if L != 0:
1205 times.append(time)
1206 energies.append([Epot, EH, Ekin])
1207 temperatures.append(temperature)
1208 pressures.append(pressure)
1209 atoms = ase.Atoms(symbols=symbols,
1210 positions=positions,
1211 cell=cell)
1212 atoms.set_velocities(velocities)
1213 if len(stress) == 0:
1214 atoms.calc = SinglePointCalculator(
1215 atoms=atoms, energy=Epot, forces=forces)
1216 else:
1217 atoms.calc = SinglePointCalculator(
1218 atoms=atoms, energy=Epot,
1219 forces=forces, stress=stress)
1220 traj.append(atoms)
1221 symbols = []
1222 positions = []
1223 velocities = []
1224 forces = []
1225 cell = []
1226 cell_velocities = []
1227 stress = []
1228 continue
1229 if len(fields) == 1:
1230 time = factors['t'] * float(fields[0])
1231 continue
1233 if fields[-1] == 'E':
1234 E = [float(x) for x in fields[0:3]]
1235 Epot, EH, Ekin = (factors['E'] * Ei for Ei in E)
1236 continue
1238 if fields[-1] == 'T':
1239 temperature = factors['T'] * float(fields[0])
1240 continue
1242 # only printed in case of variable cell calculation or calculate_stress
1243 # explicitly requested
1244 if fields[-1] == 'P':
1245 pressure = factors['P'] * float(fields[0])
1246 continue
1247 if fields[-1] == 'h':
1248 h = [float(x) for x in fields[0:3]]
1249 cell.append([factors['h'] * hi for hi in h])
1250 continue
1252 # only printed in case of variable cell calculation
1253 if fields[-1] == 'hv':
1254 hv = [float(x) for x in fields[0:3]]
1255 cell_velocities.append([factors['hv'] * hvi for hvi in hv])
1256 continue
1258 # only printed in case of variable cell calculation
1259 if fields[-1] == 'S':
1260 S = [float(x) for x in fields[0:3]]
1261 stress.append([factors['S'] * Si for Si in S])
1262 continue
1263 if fields[-1] == 'R':
1264 symbols.append(fields[0])
1265 R = [float(x) for x in fields[2:5]]
1266 positions.append([factors['R'] * Ri for Ri in R])
1267 continue
1268 if fields[-1] == 'V':
1269 V = [float(x) for x in fields[2:5]]
1270 velocities.append([factors['V'] * Vi for Vi in V])
1271 continue
1272 if fields[-1] == 'F':
1273 F = [float(x) for x in fields[2:5]]
1274 forces.append([factors['F'] * Fi for Fi in F])
1275 continue
1277 if index is None:
1278 pass
1279 else:
1280 traj = traj[index]
1282 if return_scalars:
1283 data = [times, energies, temperatures, pressures]
1284 return data, traj
1285 else:
1286 return traj
1289# Routines that only the calculator requires
1291def read_param(filename='', calc=None, fd=None, get_interface_options=False):
1292 if fd is None:
1293 if filename == '':
1294 raise ValueError('One between filename and fd must be provided')
1295 fd = open(filename)
1296 elif filename:
1297 warnings.warn('Filestream used to read param, file name will be '
1298 'ignored')
1300 # If necessary, get the interface options
1301 if get_interface_options:
1302 int_opts = {}
1303 optre = re.compile(r'# ASE_INTERFACE ([^\s]+) : ([^\s]+)')
1305 lines = fd.readlines()
1306 fd.seek(0)
1308 for line in lines:
1309 m = optre.search(line)
1310 if m:
1311 int_opts[m.groups()[0]] = m.groups()[1]
1313 data = read_freeform(fd)
1315 if calc is None:
1316 from ase.calculators.castep import Castep
1317 calc = Castep(check_castep_version=False, keyword_tolerance=2)
1319 for kw, (val, otype) in data.items():
1320 if otype == 'block':
1321 val = val.split('\n') # Avoids a bug for one-line blocks
1322 calc.param.__setattr__(kw, val)
1324 if not get_interface_options:
1325 return calc
1326 else:
1327 return calc, int_opts
1330def write_param(filename, param, check_checkfile=False,
1331 force_write=False,
1332 interface_options=None):
1333 """Writes a CastepParam object to a CASTEP .param file
1335 Parameters:
1336 filename: the location of the file to write to. If it
1337 exists it will be overwritten without warning. If it
1338 doesn't it will be created.
1339 param: a CastepParam instance
1340 check_checkfile : if set to True, write_param will
1341 only write continuation or reuse statement
1342 if a restart file exists in the same directory
1343 """
1344 if os.path.isfile(filename) and not force_write:
1345 warnings.warn('ase.io.castep.write_param: Set optional argument '
1346 'force_write=True to overwrite %s.' % filename)
1347 return False
1349 out = paropen(filename, 'w')
1350 out.write('#######################################################\n')
1351 out.write(f'#CASTEP param file: {filename}\n')
1352 out.write('#Created using the Atomic Simulation Environment (ASE)#\n')
1353 if interface_options is not None:
1354 out.write('# Internal settings of the calculator\n')
1355 out.write('# This can be switched off by settings\n')
1356 out.write('# calc._export_settings = False\n')
1357 out.write('# If stated, this will be automatically processed\n')
1358 out.write('# by ase.io.castep.read_seed()\n')
1359 for option, value in sorted(interface_options.items()):
1360 out.write(f'# ASE_INTERFACE {option} : {value}\n')
1361 out.write('#######################################################\n\n')
1363 if check_checkfile:
1364 param = deepcopy(param) # To avoid modifying the parent one
1365 for checktype in ['continuation', 'reuse']:
1366 opt = getattr(param, checktype)
1367 if opt and opt.value:
1368 fname = opt.value
1369 if fname == 'default':
1370 fname = os.path.splitext(filename)[0] + '.check'
1371 if not (os.path.exists(fname) or
1372 # CASTEP also understands relative path names, hence
1373 # also check relative to the param file directory
1374 os.path.exists(
1375 os.path.join(os.path.dirname(filename),
1376 opt.value))):
1377 opt.clear()
1379 write_freeform(out, param)
1381 out.close()
1384def read_seed(seed, new_seed=None, ignore_internal_keys=False):
1385 """A wrapper around the CASTEP Calculator in conjunction with
1386 read_cell and read_param. Basically this can be used to reuse
1387 a previous calculation which results in a triple of
1388 cell/param/castep file. The label of the calculation if pre-
1389 fixed with `copy_of_` and everything else will be recycled as
1390 much as possible from the addressed calculation.
1392 Please note that this routine will return an atoms ordering as specified
1393 in the cell file! It will thus undo the potential reordering internally
1394 done by castep.
1395 """
1397 directory = os.path.abspath(os.path.dirname(seed))
1398 seed = os.path.basename(seed)
1400 paramfile = os.path.join(directory, f'{seed}.param')
1401 cellfile = os.path.join(directory, f'{seed}.cell')
1402 castepfile = os.path.join(directory, f'{seed}.castep')
1403 checkfile = os.path.join(directory, f'{seed}.check')
1405 atoms = read_cell(cellfile)
1406 atoms.calc._directory = directory
1407 atoms.calc._rename_existing_dir = False
1408 atoms.calc._castep_pp_path = directory
1409 atoms.calc.merge_param(paramfile,
1410 ignore_internal_keys=ignore_internal_keys)
1411 if new_seed is None:
1412 atoms.calc._label = f'copy_of_{seed}'
1413 else:
1414 atoms.calc._label = str(new_seed)
1415 if os.path.isfile(castepfile):
1416 # _set_atoms needs to be True here
1417 # but we set it right back to False
1418 # atoms.calc._set_atoms = False
1419 # BUGFIX: I do not see a reason to do that!
1420 atoms.calc.read(castepfile)
1421 # atoms.calc._set_atoms = False
1423 # if here is a check file, we also want to re-use this information
1424 if os.path.isfile(checkfile):
1425 atoms.calc._check_file = os.path.basename(checkfile)
1427 # sync the top-level object with the
1428 # one attached to the calculator
1429 atoms = atoms.calc.atoms
1430 else:
1431 # There are cases where we only want to restore a calculator/atoms
1432 # setting without a castep file...
1433 # No print statement required in these cases
1434 warnings.warn(
1435 'Corresponding *.castep file not found. '
1436 'Atoms object will be restored from *.cell and *.param only.')
1437 atoms.calc.push_oldstate()
1439 return atoms
1442def read_bands(filename='', fd=None, units=units_CODATA2002):
1443 """Read Castep.bands file to kpoints, weights and eigenvalues
1445 Args:
1446 filename (str):
1447 path to seedname.bands file
1448 fd (fd):
1449 file descriptor for open bands file
1450 units (dict):
1451 Conversion factors for atomic units
1453 Returns:
1454 (tuple):
1455 (kpts, weights, eigenvalues, efermi)
1457 Where ``kpts`` and ``weights`` are 1d numpy arrays, eigenvalues
1458 is an array of shape (spin, kpts, nbands) and efermi is a float
1459 """
1461 Hartree = units['Eh']
1463 if fd is None:
1464 if filename == '':
1465 raise ValueError('One between filename and fd must be provided')
1466 fd = open(filename)
1467 elif filename:
1468 warnings.warn('Filestream used to read param, file name will be '
1469 'ignored')
1471 nkpts, nspin, _, nbands, efermi = (t(fd.readline().split()[-1]) for t in
1472 [int, int, float, int, float])
1474 kpts, weights = np.zeros((nkpts, 3)), np.zeros(nkpts)
1475 eigenvalues = np.zeros((nspin, nkpts, nbands))
1477 # Skip unit cell
1478 for _ in range(4):
1479 fd.readline()
1481 def _kptline_to_i_k_wt(line):
1482 line = line.split()
1483 line = [int(line[1])] + list(map(float, line[2:]))
1484 return (line[0] - 1, line[1:4], line[4])
1486 # CASTEP often writes these out-of-order, so check index and write directly
1487 # to the correct row
1488 for kpt_line in range(nkpts):
1489 i_kpt, kpt, wt = _kptline_to_i_k_wt(fd.readline())
1490 kpts[i_kpt, :], weights[i_kpt] = kpt, wt
1491 for spin in range(nspin):
1492 fd.readline() # Skip 'Spin component N' line
1493 eigenvalues[spin, i_kpt, :] = [float(fd.readline())
1494 for _ in range(nbands)]
1496 return (kpts, weights, eigenvalues * Hartree, efermi * Hartree)