Coverage for /builds/kinetik161/ase/ase/io/castep.py: 33.42%
751 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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
11import numpy as np
13import ase
14# independent unit management included here:
15# When high accuracy is required, this allows to easily pin down
16# unit conversion factors from different "unit definition systems"
17# (CODATA1986 for ase-3.6.0.2515 vs CODATA2002 for CASTEP 5.01).
18#
19# ase.units in in ase-3.6.0.2515 is based on CODATA1986
20import ase.units
21from ase.constraints import FixAtoms, FixCartesian, FixedLine, FixedPlane
22from ase.geometry.cell import cellpar_to_cell
23from ase.parallel import paropen
24from ase.spacegroup import Spacegroup
25from ase.utils import atoms_to_spglib_cell
27units_ase = {
28 'hbar': ase.units._hbar * ase.units.J,
29 'Eh': ase.units.Hartree,
30 'kB': ase.units.kB,
31 'a0': ase.units.Bohr,
32 't0': ase.units._hbar * ase.units.J / ase.units.Hartree,
33 'c': ase.units._c,
34 'me': ase.units._me / ase.units._amu,
35 'Pascal': 1.0 / ase.units.Pascal}
37# CODATA1986 (included herein for the sake of completeness)
38# taken from
39# http://physics.nist.gov/cuu/Archive/1986RMP.pdf
40units_CODATA1986 = {
41 'hbar': 6.5821220E-16, # eVs
42 'Eh': 27.2113961, # eV
43 'kB': 8.617385E-5, # eV/K
44 'a0': 0.529177249, # A
45 'c': 299792458, # m/s
46 'e': 1.60217733E-19, # C
47 'me': 5.485799110E-4} # u
49# CODATA2002: default in CASTEP 5.01
50# (-> check in more recent CASTEP in case of numerical discrepancies?!)
51# taken from
52# http://physics.nist.gov/cuu/Document/all_2002.pdf
53units_CODATA2002 = {
54 'hbar': 6.58211915E-16, # eVs
55 'Eh': 27.2113845, # eV
56 'kB': 8.617343E-5, # eV/K
57 'a0': 0.5291772108, # A
58 'c': 299792458, # m/s
59 'e': 1.60217653E-19, # C
60 'me': 5.4857990945E-4} # u
62# (common) derived entries
63for d in (units_CODATA1986, units_CODATA2002):
64 d['t0'] = d['hbar'] / d['Eh'] # s
65 d['Pascal'] = d['e'] * 1E30 # Pa
68__all__ = [
69 # routines for the generic io function
70 'read_castep',
71 'read_castep_castep',
72 'read_castep_castep_old',
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 ]
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))
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')
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.
126 Note that this is function is intended to maintain backwards-compatibility
127 only.
128 """
129 from ase.io import write
131 write(filename, atoms, positions_frac=positions_frac,
132 castep_cell=castep_cell, force_write=force_write)
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.
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
147 Note that the "force_write" keywords has no effect currently.
149 Arguments:
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 """
165 if atoms is None:
166 warnings.warn('Atoms object not initialized')
167 return False
168 if isinstance(atoms, list):
169 if len(atoms) > 1:
170 atoms = atoms[-1]
172 # Header
173 fd.write('#######################################################\n')
174 fd.write(f'#CASTEP cell file: {fd.name}\n')
175 fd.write('#Created using the Atomic Simulation Environment (ASE)#\n')
176 fd.write('#######################################################\n\n')
178 # To write this we simply use the existing Castep calculator, or create
179 # one
180 from ase.calculators.castep import Castep, CastepCell
182 try:
183 has_cell = isinstance(atoms.calc.cell, CastepCell)
184 except AttributeError:
185 has_cell = False
187 if has_cell:
188 cell = deepcopy(atoms.calc.cell)
189 else:
190 cell = Castep(keyword_tolerance=2).cell
192 # Write lattice
193 fformat = f'%{precision + 3}.{precision}f'
194 cell_block_format = ' '.join([fformat] * 3)
195 cell.lattice_cart = [cell_block_format % tuple(line)
196 for line in atoms.get_cell()]
198 if positions_frac:
199 pos_keyword = 'positions_frac'
200 positions = atoms.get_scaled_positions()
201 else:
202 pos_keyword = 'positions_abs'
203 positions = atoms.get_positions()
205 if atoms.has('castep_custom_species'):
206 elems = atoms.get_array('castep_custom_species')
207 else:
208 elems = atoms.get_chemical_symbols()
209 if atoms.has('masses'):
211 from ase.data import atomic_masses
212 masses = atoms.get_array('masses')
213 custom_masses = {}
215 for i, species in enumerate(elems):
216 custom_mass = masses[i]
218 # build record of different masses for each species
219 if species not in custom_masses.keys():
221 # build dictionary of positions of all species with
222 # same name and mass value ideally there should only
223 # be one mass per species
224 custom_masses[species] = {custom_mass: [i]}
226 # if multiple masses found for a species
227 elif custom_mass not in custom_masses[species].keys():
229 # if custom species were already manually defined raise an error
230 if atoms.has('castep_custom_species'):
231 raise ValueError(
232 "Could not write custom mass block for {0}. \n"
233 "Custom mass was set ({1}), but an inconsistent set of "
234 "castep_custom_species already defines "
235 "({2}) for {0}. \n"
236 "If using both features, ensure that "
237 "each species type in "
238 "atoms.arrays['castep_custom_species'] "
239 "has consistent mass values and that each atom "
240 "with non-standard "
241 "mass belongs to a custom species type."
242 "".format(
243 species, custom_mass, list(
244 custom_masses[species].keys())[0]))
246 # append mass to create custom species later
247 else:
248 custom_masses[species][custom_mass] = [i]
249 else:
250 custom_masses[species][custom_mass].append(i)
252 # create species_mass block
253 mass_block = []
255 for el, mass_dict in custom_masses.items():
257 # ignore mass record that match defaults
258 default = mass_dict.pop(atomic_masses[atoms.get_array(
259 'numbers')[list(elems).index(el)]], None)
260 if mass_dict:
261 # no custom species need to be created
262 if len(mass_dict) == 1 and not default:
263 mass_block.append('{} {}'.format(
264 el, list(mass_dict.keys())[0]))
265 # for each custom mass, create new species and change names to
266 # match in 'elems' list
267 else:
268 warnings.warn(
269 'Custom mass specified for '
270 'standard species {}, creating custom species'
271 .format(el))
273 for i, vals in enumerate(mass_dict.items()):
274 mass_val, idxs = vals
275 custom_species_name = f"{el}:{i}"
276 warnings.warn(
277 'Creating custom species {} with mass {}'.format(
278 custom_species_name, str(mass_dict)))
279 for idx in idxs:
280 elems[idx] = custom_species_name
281 mass_block.append('{} {}'.format(
282 custom_species_name, mass_val))
284 setattr(cell, 'species_mass', mass_block)
286 if atoms.has('castep_labels'):
287 labels = atoms.get_array('castep_labels')
288 else:
289 labels = ['NULL'] * len(elems)
291 if str(magnetic_moments).lower() == 'initial':
292 magmoms = atoms.get_initial_magnetic_moments()
293 elif str(magnetic_moments).lower() == 'calculated':
294 magmoms = atoms.get_magnetic_moments()
295 elif np.array(magnetic_moments).shape == (len(elems),):
296 magmoms = np.array(magnetic_moments)
297 else:
298 magmoms = [0] * len(elems)
300 pos_block = []
301 pos_block_format = '%s ' + cell_block_format
303 for i, el in enumerate(elems):
304 xyz = positions[i]
305 line = pos_block_format % tuple([el] + list(xyz))
306 # ADD other keywords if necessary
307 if magmoms[i] != 0:
308 line += f' SPIN={magmoms[i]} '
309 if labels[i].strip() not in ('NULL', ''):
310 line += f' LABEL={labels[i]} '
311 pos_block.append(line)
313 setattr(cell, pos_keyword, pos_block)
315 constraints = atoms.constraints
316 if len(constraints):
317 _supported_constraints = (FixAtoms, FixedPlane, FixedLine,
318 FixCartesian)
320 constr_block = []
322 for constr in constraints:
323 if not isinstance(constr, _supported_constraints):
324 warnings.warn(
325 'Warning: you have constraints in your atoms, that are '
326 'not supported by the CASTEP ase interface')
327 break
328 species_indices = atoms.symbols.species_indices()
329 if isinstance(constr, FixAtoms):
330 for i in constr.index:
331 try:
332 symbol = atoms.get_chemical_symbols()[i]
333 nis = species_indices[i] + 1
334 except KeyError:
335 raise RuntimeError('Unrecognized index in'
336 + f' constraint {constr}')
337 for j in range(3):
338 L = '%6d %3s %3d ' % (len(constr_block) + 1,
339 symbol,
340 nis)
341 L += ['1 0 0', '0 1 0', '0 0 1'][j]
342 constr_block += [L]
344 elif isinstance(constr, FixCartesian):
345 n = constr.a
346 symbol = atoms.get_chemical_symbols()[n]
347 nis = species_indices[n] + 1
349 for i, m in enumerate(constr.mask):
350 if m == 1:
351 continue
352 L = '%6d %3s %3d ' % (len(constr_block) + 1, symbol, nis)
353 L += ' '.join(['1' if j == i else '0' for j in range(3)])
354 constr_block += [L]
356 elif isinstance(constr, FixedPlane):
357 n = constr.a
358 symbol = atoms.get_chemical_symbols()[n]
359 nis = species_indices[n] + 1
361 L = '%6d %3s %3d ' % (len(constr_block) + 1, symbol, nis)
362 L += ' '.join([str(d) for d in constr.dir])
363 constr_block += [L]
365 elif isinstance(constr, FixedLine):
366 n = constr.a
367 symbol = atoms.get_chemical_symbols()[n]
368 nis = species_indices[n] + 1
370 direction = constr.dir
371 ((i1, v1), (i2, v2)) = sorted(enumerate(direction),
372 key=lambda x: abs(x[1]),
373 reverse=True)[:2]
374 n1 = np.zeros(3)
375 n1[i2] = v1
376 n1[i1] = -v2
377 n1 = n1 / np.linalg.norm(n1)
379 n2 = np.cross(direction, n1)
381 l1 = '%6d %3s %3d %f %f %f' % (len(constr_block) + 1,
382 symbol, nis,
383 n1[0], n1[1], n1[2])
384 l2 = '%6d %3s %3d %f %f %f' % (len(constr_block) + 2,
385 symbol, nis,
386 n2[0], n2[1], n2[2])
388 constr_block += [l1, l2]
390 cell.ionic_constraints = constr_block
392 write_freeform(fd, cell)
394 return True
397def read_freeform(fd):
398 """
399 Read a CASTEP freeform file (the basic format of .cell and .param files)
400 and return keyword-value pairs as a dict (values are strings for single
401 keywords and lists of strings for blocks).
402 """
404 from ase.calculators.castep import CastepInputFile
406 inputobj = CastepInputFile(keyword_tolerance=2)
408 filelines = fd.readlines()
410 keyw = None
411 read_block = False
412 block_lines = None
414 for i, l in enumerate(filelines):
416 # Strip all comments, aka anything after a hash
417 L = re.split(r'[#!;]', l, 1)[0].strip()
419 if L == '':
420 # Empty line... skip
421 continue
423 lsplit = re.split(r'\s*[:=]*\s+', L, 1)
425 if read_block:
426 if lsplit[0].lower() == '%endblock':
427 if len(lsplit) == 1 or lsplit[1].lower() != keyw:
428 raise ValueError('Out of place end of block at '
429 'line %i in freeform file' % i + 1)
430 else:
431 read_block = False
432 inputobj.__setattr__(keyw, block_lines)
433 else:
434 block_lines += [L]
435 else:
436 # Check the first word
438 # Is it a block?
439 read_block = (lsplit[0].lower() == '%block')
440 if read_block:
441 if len(lsplit) == 1:
442 raise ValueError(('Unrecognizable block at line %i '
443 'in io freeform file') % i + 1)
444 else:
445 keyw = lsplit[1].lower()
446 else:
447 keyw = lsplit[0].lower()
449 # Now save the value
450 if read_block:
451 block_lines = []
452 else:
453 inputobj.__setattr__(keyw, ' '.join(lsplit[1:]))
455 return inputobj.get_attr_dict(types=True)
458def read_cell(filename, index=None):
459 """
460 Wrapper function for the more generic read() functionality.
462 Note that this is function is intended to maintain backwards-compatibility
463 only.
464 """
465 from ase.io import read
466 return read(filename, index=index, format='castep-cell')
469def read_castep_cell(fd, index=None, calculator_args={}, find_spg=False,
470 units=units_CODATA2002):
471 """Read a .cell file and return an atoms object.
472 Any value found that does not fit the atoms API
473 will be stored in the atoms.calc attribute.
475 By default, the Castep calculator will be tolerant and in the absence of a
476 castep_keywords.json file it will just accept all keywords that aren't
477 automatically parsed.
478 """
480 from ase.calculators.castep import Castep
482 cell_units = { # Units specifiers for CASTEP
483 'bohr': units_CODATA2002['a0'],
484 'ang': 1.0,
485 'm': 1e10,
486 'cm': 1e8,
487 'nm': 10,
488 'pm': 1e-2
489 }
491 calc = Castep(**calculator_args)
493 if calc.cell.castep_version == 0 and calc._kw_tol < 3:
494 # No valid castep_keywords.json was found
495 warnings.warn(
496 'read_cell: Warning - Was not able to validate CASTEP input. '
497 'This may be due to a non-existing '
498 '"castep_keywords.json" '
499 'file or a non-existing CASTEP installation. '
500 'Parsing will go on but keywords will not be '
501 'validated and may cause problems if incorrect during a CASTEP '
502 'run.')
504 celldict = read_freeform(fd)
506 def parse_blockunit(line_tokens, blockname):
507 u = 1.0
508 if len(line_tokens[0]) == 1:
509 usymb = line_tokens[0][0].lower()
510 u = cell_units.get(usymb, 1)
511 if usymb not in cell_units:
512 warnings.warn('read_cell: Warning - ignoring invalid '
513 'unit specifier in %BLOCK {} '
514 '(assuming Angstrom instead)'.format(blockname))
515 line_tokens = line_tokens[1:]
516 return u, line_tokens
518 # Arguments to pass to the Atoms object at the end
519 aargs = {
520 'pbc': True
521 }
523 # Start by looking for the lattice
524 lat_keywords = [w in celldict for w in ('lattice_cart', 'lattice_abc')]
525 if all(lat_keywords):
526 warnings.warn('read_cell: Warning - two lattice blocks present in the'
527 ' same file. LATTICE_ABC will be ignored')
528 elif not any(lat_keywords):
529 raise ValueError('Cell file must contain at least one between '
530 'LATTICE_ABC and LATTICE_CART')
532 if 'lattice_abc' in celldict:
534 lines = celldict.pop('lattice_abc')[0].split('\n')
535 line_tokens = [line.split() for line in lines]
537 u, line_tokens = parse_blockunit(line_tokens, 'lattice_abc')
539 if len(line_tokens) != 2:
540 warnings.warn('read_cell: Warning - ignoring additional '
541 'lines in invalid %BLOCK LATTICE_ABC')
543 abc = [float(p) * u for p in line_tokens[0][:3]]
544 angles = [float(phi) for phi in line_tokens[1][:3]]
546 aargs['cell'] = cellpar_to_cell(abc + angles)
548 if 'lattice_cart' in celldict:
550 lines = celldict.pop('lattice_cart')[0].split('\n')
551 line_tokens = [line.split() for line in lines]
553 u, line_tokens = parse_blockunit(line_tokens, 'lattice_cart')
555 if len(line_tokens) != 3:
556 warnings.warn('read_cell: Warning - ignoring more than '
557 'three lattice vectors in invalid %BLOCK '
558 'LATTICE_CART')
560 aargs['cell'] = [[float(x) * u for x in lt[:3]] for lt in line_tokens]
562 # Now move on to the positions
563 pos_keywords = [w in celldict
564 for w in ('positions_abs', 'positions_frac')]
566 if all(pos_keywords):
567 warnings.warn('read_cell: Warning - two lattice blocks present in the'
568 ' same file. POSITIONS_FRAC will be ignored')
569 del celldict['positions_frac']
570 elif not any(pos_keywords):
571 raise ValueError('Cell file must contain at least one between '
572 'POSITIONS_FRAC and POSITIONS_ABS')
574 aargs['symbols'] = []
575 pos_type = 'positions'
576 pos_block = celldict.pop('positions_abs', [None])[0]
577 if pos_block is None:
578 pos_type = 'scaled_positions'
579 pos_block = celldict.pop('positions_frac', [None])[0]
580 aargs[pos_type] = []
582 lines = pos_block.split('\n')
583 line_tokens = [line.split() for line in lines]
585 if 'scaled' not in pos_type:
586 u, line_tokens = parse_blockunit(line_tokens, 'positions_abs')
587 else:
588 u = 1.0
590 # Here we extract all the possible additional info
591 # These are marked by their type
593 add_info = {
594 'SPIN': (float, 0.0), # (type, default)
595 'MAGMOM': (float, 0.0),
596 'LABEL': (str, 'NULL')
597 }
598 add_info_arrays = {k: [] for k in add_info}
600 def parse_info(raw_info):
602 re_keys = (r'({0})\s*[=:\s]{{1}}\s'
603 r'*([^\s]*)').format('|'.join(add_info.keys()))
604 # Capture all info groups
605 info = re.findall(re_keys, raw_info)
606 info = {g[0]: add_info[g[0]][0](g[1]) for g in info}
607 return info
609 # Array for custom species (a CASTEP special thing)
610 # Usually left unused
611 custom_species = None
613 for tokens in line_tokens:
614 # Now, process the whole 'species' thing
615 spec_custom = tokens[0].split(':', 1)
616 elem = spec_custom[0]
617 if len(spec_custom) > 1 and custom_species is None:
618 # Add it to the custom info!
619 custom_species = list(aargs['symbols'])
620 if custom_species is not None:
621 custom_species.append(tokens[0])
622 aargs['symbols'].append(elem)
623 aargs[pos_type].append([float(p) * u for p in tokens[1:4]])
624 # Now for the additional information
625 info = ' '.join(tokens[4:])
626 info = parse_info(info)
627 for k in add_info:
628 add_info_arrays[k] += [info.get(k, add_info[k][1])]
630 # read in custom species mass
631 if 'species_mass' in celldict:
632 spec_list = custom_species if custom_species else aargs['symbols']
633 aargs['masses'] = [None for _ in spec_list]
634 lines = celldict.pop('species_mass')[0].split('\n')
635 line_tokens = [line.split() for line in lines]
637 if len(line_tokens[0]) == 1:
638 if line_tokens[0][0].lower() not in ('amu', 'u'):
639 raise ValueError(
640 "unit specifier '{}' in %BLOCK SPECIES_MASS "
641 "not recognised".format(
642 line_tokens[0][0].lower()))
643 line_tokens = line_tokens[1:]
645 for tokens in line_tokens:
646 token_pos_list = [i for i, x in enumerate(
647 spec_list) if x == tokens[0]]
648 if len(token_pos_list) == 0:
649 warnings.warn(
650 'read_cell: Warning - ignoring unused '
651 'species mass {} in %BLOCK SPECIES_MASS'.format(
652 tokens[0]))
653 for idx in token_pos_list:
654 aargs['masses'][idx] = tokens[1]
656 # Now on to the species potentials...
657 if 'species_pot' in celldict:
658 lines = celldict.pop('species_pot')[0].split('\n')
659 line_tokens = [line.split() for line in lines]
661 for tokens in line_tokens:
662 if len(tokens) == 1:
663 # It's a library
664 all_spec = (set(custom_species) if custom_species is not None
665 else set(aargs['symbols']))
666 for s in all_spec:
667 calc.cell.species_pot = (s, tokens[0])
668 else:
669 calc.cell.species_pot = tuple(tokens[:2])
671 # Ionic constraints
672 raw_constraints = {}
674 if 'ionic_constraints' in celldict:
675 lines = celldict.pop('ionic_constraints')[0].split('\n')
676 line_tokens = [line.split() for line in lines]
678 for tokens in line_tokens:
679 if not len(tokens) == 6:
680 continue
681 _, species, nic, x, y, z = tokens
682 # convert xyz to floats
683 x = float(x)
684 y = float(y)
685 z = float(z)
687 nic = int(nic)
688 if (species, nic) not in raw_constraints:
689 raw_constraints[(species, nic)] = []
690 raw_constraints[(species, nic)].append(np.array(
691 [x, y, z]))
693 # Symmetry operations
694 if 'symmetry_ops' in celldict:
695 lines = celldict.pop('symmetry_ops')[0].split('\n')
696 line_tokens = [line.split() for line in lines]
698 # Read them in blocks of four
699 blocks = np.array(line_tokens).astype(float)
700 if (len(blocks.shape) != 2 or blocks.shape[1] != 3
701 or blocks.shape[0] % 4 != 0):
702 warnings.warn('Warning: could not parse SYMMETRY_OPS'
703 ' block properly, skipping')
704 else:
705 blocks = blocks.reshape((-1, 4, 3))
706 rotations = blocks[:, :3]
707 translations = blocks[:, 3]
709 # Regardless of whether we recognize them, store these
710 calc.cell.symmetry_ops = (rotations, translations)
712 # Anything else that remains, just add it to the cell object:
713 for k, (val, otype) in celldict.items():
714 try:
715 if otype == 'block':
716 val = val.split('\n') # Avoids a bug for one-line blocks
717 calc.cell.__setattr__(k, val)
718 except Exception as e:
719 raise RuntimeError(
720 f'Problem setting calc.cell.{k} = {val}: {e}')
722 # Get the relevant additional info
723 aargs['magmoms'] = np.array(add_info_arrays['SPIN'])
724 # SPIN or MAGMOM are alternative keywords
725 aargs['magmoms'] = np.where(aargs['magmoms'] != 0,
726 aargs['magmoms'],
727 add_info_arrays['MAGMOM'])
728 labels = np.array(add_info_arrays['LABEL'])
730 aargs['calculator'] = calc
732 atoms = ase.Atoms(**aargs)
734 # Spacegroup...
735 if find_spg:
736 # Try importing spglib
737 try:
738 import spglib
739 except ImportError:
740 warnings.warn('spglib not found installed on this system - '
741 'automatic spacegroup detection is not possible')
742 spglib = None
744 if spglib is not None:
745 symmd = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms))
746 atoms_spg = Spacegroup(int(symmd['number']))
747 atoms.info['spacegroup'] = atoms_spg
749 atoms.new_array('castep_labels', labels)
750 if custom_species is not None:
751 atoms.new_array('castep_custom_species', np.array(custom_species))
753 fixed_atoms = []
754 constraints = []
755 index_dict = atoms.symbols.indices()
756 for (species, nic), value in raw_constraints.items():
758 absolute_nr = index_dict[species][nic - 1]
759 if len(value) == 3:
760 # Check if they are linearly independent
761 if np.linalg.det(value) == 0:
762 warnings.warn(
763 'Error: Found linearly dependent constraints attached '
764 'to atoms %s' %
765 (absolute_nr))
766 continue
767 fixed_atoms.append(absolute_nr)
768 elif len(value) == 2:
769 direction = np.cross(value[0], value[1])
770 # Check if they are linearly independent
771 if np.linalg.norm(direction) == 0:
772 warnings.warn(
773 'Error: Found linearly dependent constraints attached '
774 'to atoms %s' %
775 (absolute_nr))
776 continue
777 constraint = ase.constraints.FixedLine(
778 a=absolute_nr,
779 direction=direction)
780 constraints.append(constraint)
781 elif len(value) == 1:
782 constraint = ase.constraints.FixedPlane(
783 a=absolute_nr,
784 direction=np.array(value[0], dtype=np.float32))
785 constraints.append(constraint)
786 else:
787 warnings.warn(
788 f'Error: Found {len(value)} statements attached to atoms '
789 f'{absolute_nr}'
790 )
792 # we need to sort the fixed atoms list in order not to raise an assertion
793 # error in FixAtoms
794 if fixed_atoms:
795 constraints.append(
796 ase.constraints.FixAtoms(indices=sorted(fixed_atoms)))
797 if constraints:
798 atoms.set_constraint(constraints)
800 atoms.calc.atoms = atoms
801 atoms.calc.push_oldstate()
803 return atoms
806def read_castep(filename, index=None):
807 """
808 Wrapper function for the more generic read() functionality.
810 Note that this is function is intended to maintain backwards-compatibility
811 only.
812 """
813 from ase.io import read
814 return read(filename, index=index, format='castep-castep')
817def read_castep_castep(fd, index=None):
818 """
819 Reads a .castep file and returns an atoms object.
820 The calculator information will be stored in the calc attribute.
822 There is no use of the "index" argument as of now, it is just inserted for
823 convenience to comply with the generic "read()" in ase.io
825 Please note that this routine will return an atom ordering as found
826 within the castep file. This means that the species will be ordered by
827 ascending atomic numbers. The atoms witin a species are ordered as given
828 in the original cell file.
830 Note: This routine returns a single atoms_object only, the last
831 configuration in the file. Yet, if you want to parse an MD run, use the
832 novel function `read_md()`
833 """
835 from ase.calculators.castep import Castep
837 try:
838 calc = Castep()
839 except Exception as e:
840 # No CASTEP keywords found?
841 warnings.warn(f'WARNING: {e} Using fallback .castep reader...')
842 # Fall back on the old method
843 return read_castep_castep_old(fd, index)
845 calc.read(castep_file=fd)
847 # now we trick the calculator instance such that we can savely extract
848 # energies and forces from this atom. Basically what we do is to trick the
849 # internal routine calculation_required() to always return False such that
850 # we do not need to re-run a CASTEP calculation.
851 #
852 # Probably we can solve this with a flag to the read() routine at some
853 # point, but for the moment I do not want to change too much in there.
854 calc._old_atoms = calc.atoms
855 calc._old_param = calc.param
856 calc._old_cell = calc.cell
858 return [calc.atoms] # Returning in the form of a list for next()
861def read_castep_castep_old(fd, index=None):
862 """
863 DEPRECATED
864 Now replaced by ase.calculators.castep.Castep.read(). Left in for future
865 reference and backwards compatibility needs, as well as a fallback for
866 when castep_keywords.py can't be created.
868 Reads a .castep file and returns an atoms object.
869 The calculator information will be stored in the calc attribute.
870 If more than one SCF step is found, a list of all steps
871 will be stored in the traj attribute.
873 Note that the index argument has no effect as of now.
875 Please note that this routine will return an atom ordering as found
876 within the castep file. This means that the species will be ordered by
877 ascending atomic numbers. The atoms witin a species are ordered as given
878 in the original cell file.
879 """
880 from ase.calculators.singlepoint import SinglePointCalculator
882 lines = fd.readlines()
884 traj = []
885 energy_total = None
886 energy_0K = None
887 for i, line in enumerate(lines):
888 if 'NB est. 0K energy' in line:
889 energy_0K = float(line.split()[6])
890 # support also for dispersion correction
891 elif 'NB dispersion corrected est. 0K energy*' in line:
892 energy_0K = float(line.split()[-2])
893 elif 'Final energy, E' in line:
894 energy_total = float(line.split()[4])
895 elif 'Dispersion corrected final energy' in line:
896 pass
897 # dispcorr_energy_total = float(line.split()[-2])
898 # sedc_apply = True
899 elif 'Dispersion corrected final free energy' in line:
900 pass # dispcorr_energy_free = float(line.split()[-2])
901 elif 'dispersion corrected est. 0K energy' in line:
902 pass # dispcorr_energy_0K = float(line.split()[-2])
903 elif 'Unit Cell' in line:
904 cell = [x.split()[0:3] for x in lines[i + 3:i + 6]]
905 cell = np.array([[float(col) for col in row] for row in cell])
906 elif 'Cell Contents' in line:
907 geom_starts = i
908 start_found = False
909 for j, jline in enumerate(lines[geom_starts:]):
910 if jline.find('xxxxx') > 0 and start_found:
911 geom_stop = j + geom_starts
912 break
913 if jline.find('xxxx') > 0 and not start_found:
914 geom_start = j + geom_starts + 4
915 start_found = True
916 species = [line.split()[1] for line in lines[geom_start:geom_stop]]
917 geom = np.dot(np.array([[float(col) for col in line.split()[3:6]]
918 for line in lines[geom_start:geom_stop]]),
919 cell)
920 elif 'Writing model to' in line:
921 atoms = ase.Atoms(
922 cell=cell,
923 pbc=True,
924 positions=geom,
925 symbols=''.join(species))
926 # take 0K energy where available, else total energy
927 if energy_0K:
928 energy = energy_0K
929 else:
930 energy = energy_total
931 # generate a minimal single-point calculator
932 sp_calc = SinglePointCalculator(atoms=atoms,
933 energy=energy,
934 forces=None,
935 magmoms=None,
936 stress=None)
937 atoms.calc = sp_calc
938 traj.append(atoms)
939 if index is None:
940 return traj
941 else:
942 return traj[index]
945def read_geom(filename, index=':', units=units_CODATA2002):
946 """
947 Wrapper function for the more generic read() functionality.
949 Note that this is function is intended to maintain backwards-compatibility
950 only. Keyword arguments will be passed to read_castep_geom().
951 """
952 from ase.io import read
953 return read(filename, index=index, format='castep-geom', units=units)
956def read_castep_geom(fd, index=None, units=units_CODATA2002):
957 """Reads a .geom file produced by the CASTEP GeometryOptimization task and
958 returns an atoms object.
959 The information about total free energy and forces of each atom for every
960 relaxation step will be stored for further analysis especially in a
961 single-point calculator.
962 Note that everything in the .geom file is in atomic units, which has
963 been conversed to commonly used unit angstrom(length) and eV (energy).
965 Note that the index argument has no effect as of now.
967 Contribution by Wei-Bing Zhang. Thanks!
969 Routine now accepts a filedescriptor in order to out-source the gz and
970 bz2 handling to formats.py. Note that there is a fallback routine
971 read_geom() that behaves like previous versions did.
972 """
973 from ase.calculators.singlepoint import SinglePointCalculator
975 # fd is closed by embracing read() routine
976 txt = fd.readlines()
978 traj = []
980 Hartree = units['Eh']
981 Bohr = units['a0']
983 # Yeah, we know that...
984 # print('N.B.: Energy in .geom file is not 0K extrapolated.')
985 for i, line in enumerate(txt):
986 if line.find('<-- E') > 0:
987 start_found = True
988 energy = float(line.split()[0]) * Hartree
989 cell = [x.split()[0:3] for x in txt[i + 1:i + 4]]
990 cell = np.array([[float(col) * Bohr for col in row] for row in
991 cell])
992 if line.find('<-- R') > 0 and start_found:
993 start_found = False
994 geom_start = i
995 for i, line in enumerate(txt[geom_start:]):
996 if line.find('<-- F') > 0:
997 geom_stop = i + geom_start
998 break
999 species = [line.split()[0] for line in
1000 txt[geom_start:geom_stop]]
1001 geom = np.array([[float(col) * Bohr for col in
1002 line.split()[2:5]] for line in
1003 txt[geom_start:geom_stop]])
1004 forces = np.array([[float(col) * Hartree / Bohr for col in
1005 line.split()[2:5]] for line in
1006 txt[geom_stop:geom_stop
1007 + (geom_stop - geom_start)]])
1008 image = ase.Atoms(species, geom, cell=cell, pbc=True)
1009 image.calc = SinglePointCalculator(
1010 atoms=image, energy=energy, forces=forces)
1011 traj.append(image)
1013 if index is None:
1014 return traj
1015 else:
1016 return traj[index]
1019def read_phonon(filename, index=None, read_vib_data=False,
1020 gamma_only=True, frequency_factor=None,
1021 units=units_CODATA2002):
1022 """
1023 Wrapper function for the more generic read() functionality.
1025 Note that this is function is intended to maintain backwards-compatibility
1026 only. For documentation see read_castep_phonon().
1027 """
1028 from ase.io import read
1030 if read_vib_data:
1031 full_output = True
1032 else:
1033 full_output = False
1035 return read(filename, index=index, format='castep-phonon',
1036 full_output=full_output, read_vib_data=read_vib_data,
1037 gamma_only=gamma_only, frequency_factor=frequency_factor,
1038 units=units)
1041def read_castep_phonon(fd, index=None, read_vib_data=False,
1042 gamma_only=True, frequency_factor=None,
1043 units=units_CODATA2002):
1044 """
1045 Reads a .phonon file written by a CASTEP Phonon task and returns an atoms
1046 object, as well as the calculated vibrational data if requested.
1048 Note that the index argument has no effect as of now.
1049 """
1051 # fd is closed by embracing read() routine
1052 lines = fd.readlines()
1054 atoms = None
1055 cell = []
1056 N = Nb = Nq = 0
1057 scaled_positions = []
1058 symbols = []
1059 masses = []
1061 # header
1062 L = 0
1063 while L < len(lines):
1065 line = lines[L]
1067 if 'Number of ions' in line:
1068 N = int(line.split()[3])
1069 elif 'Number of branches' in line:
1070 Nb = int(line.split()[3])
1071 elif 'Number of wavevectors' in line:
1072 Nq = int(line.split()[3])
1073 elif 'Unit cell vectors (A)' in line:
1074 for ll in range(3):
1075 L += 1
1076 fields = lines[L].split()
1077 cell.append([float(x) for x in fields[0:3]])
1078 elif 'Fractional Co-ordinates' in line:
1079 for ll in range(N):
1080 L += 1
1081 fields = lines[L].split()
1082 scaled_positions.append([float(x) for x in fields[1:4]])
1083 symbols.append(fields[4])
1084 masses.append(float(fields[5]))
1085 elif 'END header' in line:
1086 L += 1
1087 atoms = ase.Atoms(symbols=symbols,
1088 scaled_positions=scaled_positions,
1089 cell=cell)
1090 break
1092 L += 1
1094 # Eigenmodes and -vectors
1095 if frequency_factor is None:
1096 Kayser_to_eV = 1E2 * 2 * np.pi * units['hbar'] * units['c']
1097 # N.B. "fixed default" unit for frequencies in .phonon files is "cm-1"
1098 # (i.e. the latter is unaffected by the internal unit conversion system of
1099 # CASTEP!) set conversion factor to convert therefrom to eV by default for
1100 # now
1101 frequency_factor = Kayser_to_eV
1102 qpoints = []
1103 weights = []
1104 frequencies = []
1105 displacements = []
1106 for nq in range(Nq):
1107 fields = lines[L].split()
1108 qpoints.append([float(x) for x in fields[2:5]])
1109 weights.append(float(fields[5]))
1110 freqs = []
1111 for ll in range(Nb):
1112 L += 1
1113 fields = lines[L].split()
1114 freqs.append(frequency_factor * float(fields[1]))
1115 frequencies.append(np.array(freqs))
1117 # skip the two Phonon Eigenvectors header lines
1118 L += 2
1120 # generate a list of displacements with a structure that is identical to
1121 # what is stored internally in the Vibrations class (see in
1122 # ase.vibrations.Vibrations.modes):
1123 # np.array(displacements).shape == (Nb,3*N)
1125 disps = []
1126 for ll in range(Nb):
1127 disp_coords = []
1128 for lll in range(N):
1129 L += 1
1130 fields = lines[L].split()
1131 disp_x = float(fields[2]) + float(fields[3]) * 1.0j
1132 disp_y = float(fields[4]) + float(fields[5]) * 1.0j
1133 disp_z = float(fields[6]) + float(fields[7]) * 1.0j
1134 disp_coords.extend([disp_x, disp_y, disp_z])
1135 disps.append(np.array(disp_coords))
1136 displacements.append(np.array(disps))
1138 if read_vib_data:
1139 if gamma_only:
1140 vibdata = [frequencies[0], displacements[0]]
1141 else:
1142 vibdata = [qpoints, weights, frequencies, displacements]
1143 return vibdata, atoms
1144 else:
1145 return atoms
1148def read_md(filename, index=None, return_scalars=False,
1149 units=units_CODATA2002):
1150 """Wrapper function for the more generic read() functionality.
1152 Note that this function is intended to maintain backwards-compatibility
1153 only. For documentation see read_castep_md()
1154 """
1155 if return_scalars:
1156 full_output = True
1157 else:
1158 full_output = False
1160 from ase.io import read
1161 return read(filename, index=index, format='castep-md',
1162 full_output=full_output, return_scalars=return_scalars,
1163 units=units)
1166def read_castep_md(fd, index=None, return_scalars=False,
1167 units=units_CODATA2002):
1168 """Reads a .md file written by a CASTEP MolecularDynamics task
1169 and returns the trajectory stored therein as a list of atoms object.
1171 Note that the index argument has no effect as of now."""
1173 from ase.calculators.singlepoint import SinglePointCalculator
1175 factors = {
1176 't': units['t0'] * 1E15, # fs
1177 'E': units['Eh'], # eV
1178 'T': units['Eh'] / units['kB'],
1179 'P': units['Eh'] / units['a0']**3 * units['Pascal'],
1180 'h': units['a0'],
1181 'hv': units['a0'] / units['t0'],
1182 'S': units['Eh'] / units['a0']**3,
1183 'R': units['a0'],
1184 'V': np.sqrt(units['Eh'] / units['me']),
1185 'F': units['Eh'] / units['a0']}
1187 # fd is closed by embracing read() routine
1188 lines = fd.readlines()
1190 L = 0
1191 while 'END header' not in lines[L]:
1192 L += 1
1193 l_end_header = L
1194 lines = lines[l_end_header + 1:]
1195 times = []
1196 energies = []
1197 temperatures = []
1198 pressures = []
1199 traj = []
1201 # Initialization
1202 time = None
1203 Epot = None
1204 Ekin = None
1205 EH = None
1206 temperature = None
1207 pressure = None
1208 symbols = None
1209 positions = None
1210 cell = None
1211 velocities = None
1212 symbols = []
1213 positions = []
1214 velocities = []
1215 forces = []
1216 cell = np.eye(3)
1217 cell_velocities = []
1218 stress = []
1220 for (L, line) in enumerate(lines):
1221 fields = line.split()
1222 if len(fields) == 0:
1223 if L != 0:
1224 times.append(time)
1225 energies.append([Epot, EH, Ekin])
1226 temperatures.append(temperature)
1227 pressures.append(pressure)
1228 atoms = ase.Atoms(symbols=symbols,
1229 positions=positions,
1230 cell=cell)
1231 atoms.set_velocities(velocities)
1232 if len(stress) == 0:
1233 atoms.calc = SinglePointCalculator(
1234 atoms=atoms, energy=Epot, forces=forces)
1235 else:
1236 atoms.calc = SinglePointCalculator(
1237 atoms=atoms, energy=Epot,
1238 forces=forces, stress=stress)
1239 traj.append(atoms)
1240 symbols = []
1241 positions = []
1242 velocities = []
1243 forces = []
1244 cell = []
1245 cell_velocities = []
1246 stress = []
1247 continue
1248 if len(fields) == 1:
1249 time = factors['t'] * float(fields[0])
1250 continue
1252 if fields[-1] == 'E':
1253 E = [float(x) for x in fields[0:3]]
1254 Epot, EH, Ekin = (factors['E'] * Ei for Ei in E)
1255 continue
1257 if fields[-1] == 'T':
1258 temperature = factors['T'] * float(fields[0])
1259 continue
1261 # only printed in case of variable cell calculation or calculate_stress
1262 # explicitly requested
1263 if fields[-1] == 'P':
1264 pressure = factors['P'] * float(fields[0])
1265 continue
1266 if fields[-1] == 'h':
1267 h = [float(x) for x in fields[0:3]]
1268 cell.append([factors['h'] * hi for hi in h])
1269 continue
1271 # only printed in case of variable cell calculation
1272 if fields[-1] == 'hv':
1273 hv = [float(x) for x in fields[0:3]]
1274 cell_velocities.append([factors['hv'] * hvi for hvi in hv])
1275 continue
1277 # only printed in case of variable cell calculation
1278 if fields[-1] == 'S':
1279 S = [float(x) for x in fields[0:3]]
1280 stress.append([factors['S'] * Si for Si in S])
1281 continue
1282 if fields[-1] == 'R':
1283 symbols.append(fields[0])
1284 R = [float(x) for x in fields[2:5]]
1285 positions.append([factors['R'] * Ri for Ri in R])
1286 continue
1287 if fields[-1] == 'V':
1288 V = [float(x) for x in fields[2:5]]
1289 velocities.append([factors['V'] * Vi for Vi in V])
1290 continue
1291 if fields[-1] == 'F':
1292 F = [float(x) for x in fields[2:5]]
1293 forces.append([factors['F'] * Fi for Fi in F])
1294 continue
1296 if index is None:
1297 pass
1298 else:
1299 traj = traj[index]
1301 if return_scalars:
1302 data = [times, energies, temperatures, pressures]
1303 return data, traj
1304 else:
1305 return traj
1308# Routines that only the calculator requires
1310def read_param(filename='', calc=None, fd=None, get_interface_options=False):
1311 if fd is None:
1312 if filename == '':
1313 raise ValueError('One between filename and fd must be provided')
1314 fd = open(filename)
1315 elif filename:
1316 warnings.warn('Filestream used to read param, file name will be '
1317 'ignored')
1319 # If necessary, get the interface options
1320 if get_interface_options:
1321 int_opts = {}
1322 optre = re.compile(r'# ASE_INTERFACE ([^\s]+) : ([^\s]+)')
1324 lines = fd.readlines()
1325 fd.seek(0)
1327 for line in lines:
1328 m = optre.search(line)
1329 if m:
1330 int_opts[m.groups()[0]] = m.groups()[1]
1332 data = read_freeform(fd)
1334 if calc is None:
1335 from ase.calculators.castep import Castep
1336 calc = Castep(check_castep_version=False, keyword_tolerance=2)
1338 for kw, (val, otype) in data.items():
1339 if otype == 'block':
1340 val = val.split('\n') # Avoids a bug for one-line blocks
1341 calc.param.__setattr__(kw, val)
1343 if not get_interface_options:
1344 return calc
1345 else:
1346 return calc, int_opts
1349def write_param(filename, param, check_checkfile=False,
1350 force_write=False,
1351 interface_options=None):
1352 """Writes a CastepParam object to a CASTEP .param file
1354 Parameters:
1355 filename: the location of the file to write to. If it
1356 exists it will be overwritten without warning. If it
1357 doesn't it will be created.
1358 param: a CastepParam instance
1359 check_checkfile : if set to True, write_param will
1360 only write continuation or reuse statement
1361 if a restart file exists in the same directory
1362 """
1363 if os.path.isfile(filename) and not force_write:
1364 warnings.warn('ase.io.castep.write_param: Set optional argument '
1365 'force_write=True to overwrite %s.' % filename)
1366 return False
1368 out = paropen(filename, 'w')
1369 out.write('#######################################################\n')
1370 out.write(f'#CASTEP param file: {filename}\n')
1371 out.write('#Created using the Atomic Simulation Environment (ASE)#\n')
1372 if interface_options is not None:
1373 out.write('# Internal settings of the calculator\n')
1374 out.write('# This can be switched off by settings\n')
1375 out.write('# calc._export_settings = False\n')
1376 out.write('# If stated, this will be automatically processed\n')
1377 out.write('# by ase.io.castep.read_seed()\n')
1378 for option, value in sorted(interface_options.items()):
1379 out.write(f'# ASE_INTERFACE {option} : {value}\n')
1380 out.write('#######################################################\n\n')
1382 if check_checkfile:
1383 param = deepcopy(param) # To avoid modifying the parent one
1384 for checktype in ['continuation', 'reuse']:
1385 opt = getattr(param, checktype)
1386 if opt and opt.value:
1387 fname = opt.value
1388 if fname == 'default':
1389 fname = os.path.splitext(filename)[0] + '.check'
1390 if not (os.path.exists(fname) or
1391 # CASTEP also understands relative path names, hence
1392 # also check relative to the param file directory
1393 os.path.exists(
1394 os.path.join(os.path.dirname(filename),
1395 opt.value))):
1396 opt.clear()
1398 write_freeform(out, param)
1400 out.close()
1403def read_seed(seed, new_seed=None, ignore_internal_keys=False):
1404 """A wrapper around the CASTEP Calculator in conjunction with
1405 read_cell and read_param. Basically this can be used to reuse
1406 a previous calculation which results in a triple of
1407 cell/param/castep file. The label of the calculation if pre-
1408 fixed with `copy_of_` and everything else will be recycled as
1409 much as possible from the addressed calculation.
1411 Please note that this routine will return an atoms ordering as specified
1412 in the cell file! It will thus undo the potential reordering internally
1413 done by castep.
1414 """
1416 directory = os.path.abspath(os.path.dirname(seed))
1417 seed = os.path.basename(seed)
1419 paramfile = os.path.join(directory, f'{seed}.param')
1420 cellfile = os.path.join(directory, f'{seed}.cell')
1421 castepfile = os.path.join(directory, f'{seed}.castep')
1422 checkfile = os.path.join(directory, f'{seed}.check')
1424 atoms = read_cell(cellfile)
1425 atoms.calc._directory = directory
1426 atoms.calc._rename_existing_dir = False
1427 atoms.calc._castep_pp_path = directory
1428 atoms.calc.merge_param(paramfile,
1429 ignore_internal_keys=ignore_internal_keys)
1430 if new_seed is None:
1431 atoms.calc._label = f'copy_of_{seed}'
1432 else:
1433 atoms.calc._label = str(new_seed)
1434 if os.path.isfile(castepfile):
1435 # _set_atoms needs to be True here
1436 # but we set it right back to False
1437 # atoms.calc._set_atoms = False
1438 # BUGFIX: I do not see a reason to do that!
1439 atoms.calc.read(castepfile)
1440 # atoms.calc._set_atoms = False
1442 # if here is a check file, we also want to re-use this information
1443 if os.path.isfile(checkfile):
1444 atoms.calc._check_file = os.path.basename(checkfile)
1446 # sync the top-level object with the
1447 # one attached to the calculator
1448 atoms = atoms.calc.atoms
1449 else:
1450 # There are cases where we only want to restore a calculator/atoms
1451 # setting without a castep file...
1452 pass
1453 # No print statement required in these cases
1454 warnings.warn(
1455 'Corresponding *.castep file not found. '
1456 'Atoms object will be restored from *.cell and *.param only.')
1457 atoms.calc.push_oldstate()
1459 return atoms
1462def read_bands(filename='', fd=None, units=units_CODATA2002):
1463 """Read Castep.bands file to kpoints, weights and eigenvalues
1465 Args:
1466 filename (str):
1467 path to seedname.bands file
1468 fd (fd):
1469 file descriptor for open bands file
1470 units (dict):
1471 Conversion factors for atomic units
1473 Returns:
1474 (tuple):
1475 (kpts, weights, eigenvalues, efermi)
1477 Where ``kpts`` and ``weights`` are 1d numpy arrays, eigenvalues
1478 is an array of shape (spin, kpts, nbands) and efermi is a float
1479 """
1481 Hartree = units['Eh']
1483 if fd is None:
1484 if filename == '':
1485 raise ValueError('One between filename and fd must be provided')
1486 fd = open(filename)
1487 elif filename:
1488 warnings.warn('Filestream used to read param, file name will be '
1489 'ignored')
1491 nkpts, nspin, _, nbands, efermi = (t(fd.readline().split()[-1]) for t in
1492 [int, int, float, int, float])
1494 kpts, weights = np.zeros((nkpts, 3)), np.zeros(nkpts)
1495 eigenvalues = np.zeros((nspin, nkpts, nbands))
1497 # Skip unit cell
1498 for _ in range(4):
1499 fd.readline()
1501 def _kptline_to_i_k_wt(line):
1502 line = line.split()
1503 line = [int(line[1])] + list(map(float, line[2:]))
1504 return (line[0] - 1, line[1:4], line[4])
1506 # CASTEP often writes these out-of-order, so check index and write directly
1507 # to the correct row
1508 for kpt_line in range(nkpts):
1509 i_kpt, kpt, wt = _kptline_to_i_k_wt(fd.readline())
1510 kpts[i_kpt, :], weights[i_kpt] = kpt, wt
1511 for spin in range(nspin):
1512 fd.readline() # Skip 'Spin component N' line
1513 eigenvalues[spin, i_kpt, :] = [float(fd.readline())
1514 for _ in range(nbands)]
1516 return (kpts, weights, eigenvalues * Hartree, efermi * Hartree)