Coverage for /builds/ase/ase/ase/io/castep.py : 33.47%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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
9import numpy as np
10from copy import deepcopy
12import ase
14from ase.parallel import paropen
15from ase.spacegroup import Spacegroup
16from ase.geometry.cell import cellpar_to_cell
17from ase.constraints import FixAtoms, FixedPlane, FixedLine, FixCartesian
18from ase.utils import atoms_to_spglib_cell
20# independent unit management included here:
21# When high accuracy is required, this allows to easily pin down
22# unit conversion factors from different "unit definition systems"
23# (CODATA1986 for ase-3.6.0.2515 vs CODATA2002 for CASTEP 5.01).
24#
25# ase.units in in ase-3.6.0.2515 is based on CODATA1986
26import ase.units
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('{0}: {1}\n'.format(kw.upper(), opt.value))
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('#CASTEP cell file: %s\n' % fd.name)
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 = '%{0}.{1}f'.format(precision + 3, precision)
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('{0} {1}'.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 {0}, creating custom species'
271 .format(el))
273 for i, vals in enumerate(mass_dict.items()):
274 mass_val, idxs = vals
275 custom_species_name = "{0}:{1}".format(el, i)
276 warnings.warn(
277 'Creating custom species {0} with mass {1}'.format(
278 custom_species_name, str(mass_dict)))
279 for idx in idxs:
280 elems[idx] = custom_species_name
281 mass_block.append('{0} {1}'.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 += ' SPIN={0} '.format(magmoms[i])
309 if labels[i].strip() not in ('NULL', ''):
310 line += ' LABEL={0} '.format(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 UserWarning('Unrecognized index in'
336 + ' constraint %s' % 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 {0} '
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 = dict((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 '{0}' 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 {0} 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 'Problem setting calc.cell.%s = %s: %s' % (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('Error: Found %s statements attached to atoms %s' %
788 (len(value), absolute_nr))
790 # we need to sort the fixed atoms list in order not to raise an assertion
791 # error in FixAtoms
792 if fixed_atoms:
793 constraints.append(
794 ase.constraints.FixAtoms(indices=sorted(fixed_atoms)))
795 if constraints:
796 atoms.set_constraint(constraints)
798 atoms.calc.atoms = atoms
799 atoms.calc.push_oldstate()
801 return atoms
804def read_castep(filename, index=None):
805 """
806 Wrapper function for the more generic read() functionality.
808 Note that this is function is intended to maintain backwards-compatibility
809 only.
810 """
811 from ase.io import read
812 return read(filename, index=index, format='castep-castep')
815def read_castep_castep(fd, index=None):
816 """
817 Reads a .castep file and returns an atoms object.
818 The calculator information will be stored in the calc attribute.
820 There is no use of the "index" argument as of now, it is just inserted for
821 convenience to comply with the generic "read()" in ase.io
823 Please note that this routine will return an atom ordering as found
824 within the castep file. This means that the species will be ordered by
825 ascending atomic numbers. The atoms witin a species are ordered as given
826 in the original cell file.
828 Note: This routine returns a single atoms_object only, the last
829 configuration in the file. Yet, if you want to parse an MD run, use the
830 novel function `read_md()`
831 """
833 from ase.calculators.castep import Castep
835 try:
836 calc = Castep()
837 except Exception as e:
838 # No CASTEP keywords found?
839 warnings.warn('WARNING: {0} Using fallback .castep reader...'.format(e))
840 # Fall back on the old method
841 return read_castep_castep_old(fd, index)
843 calc.read(castep_file=fd)
845 # now we trick the calculator instance such that we can savely extract
846 # energies and forces from this atom. Basically what we do is to trick the
847 # internal routine calculation_required() to always return False such that
848 # we do not need to re-run a CASTEP calculation.
849 #
850 # Probably we can solve this with a flag to the read() routine at some
851 # point, but for the moment I do not want to change too much in there.
852 calc._old_atoms = calc.atoms
853 calc._old_param = calc.param
854 calc._old_cell = calc.cell
856 return [calc.atoms] # Returning in the form of a list for next()
859def read_castep_castep_old(fd, index=None):
860 """
861 DEPRECATED
862 Now replaced by ase.calculators.castep.Castep.read(). Left in for future
863 reference and backwards compatibility needs, as well as a fallback for
864 when castep_keywords.py can't be created.
866 Reads a .castep file and returns an atoms object.
867 The calculator information will be stored in the calc attribute.
868 If more than one SCF step is found, a list of all steps
869 will be stored in the traj attribute.
871 Note that the index argument has no effect as of now.
873 Please note that this routine will return an atom ordering as found
874 within the castep file. This means that the species will be ordered by
875 ascending atomic numbers. The atoms witin a species are ordered as given
876 in the original cell file.
877 """
878 from ase.calculators.singlepoint import SinglePointCalculator
880 lines = fd.readlines()
882 traj = []
883 energy_total = None
884 energy_0K = None
885 for i, line in enumerate(lines):
886 if 'NB est. 0K energy' in line:
887 energy_0K = float(line.split()[6])
888 # support also for dispersion correction
889 elif 'NB dispersion corrected est. 0K energy*' in line:
890 energy_0K = float(line.split()[-2])
891 elif 'Final energy, E' in line:
892 energy_total = float(line.split()[4])
893 elif 'Dispersion corrected final energy' in line:
894 pass
895 # dispcorr_energy_total = float(line.split()[-2])
896 # sedc_apply = True
897 elif 'Dispersion corrected final free energy' in line:
898 pass # dispcorr_energy_free = float(line.split()[-2])
899 elif 'dispersion corrected est. 0K energy' in line:
900 pass # dispcorr_energy_0K = float(line.split()[-2])
901 elif 'Unit Cell' in line:
902 cell = [x.split()[0:3] for x in lines[i + 3:i + 6]]
903 cell = np.array([[float(col) for col in row] for row in cell])
904 elif 'Cell Contents' in line:
905 geom_starts = i
906 start_found = False
907 for j, jline in enumerate(lines[geom_starts:]):
908 if jline.find('xxxxx') > 0 and start_found:
909 geom_stop = j + geom_starts
910 break
911 if jline.find('xxxx') > 0 and not start_found:
912 geom_start = j + geom_starts + 4
913 start_found = True
914 species = [line.split()[1] for line in lines[geom_start:geom_stop]]
915 geom = np.dot(np.array([[float(col) for col in line.split()[3:6]]
916 for line in lines[geom_start:geom_stop]]),
917 cell)
918 elif 'Writing model to' in line:
919 atoms = ase.Atoms(
920 cell=cell,
921 pbc=True,
922 positions=geom,
923 symbols=''.join(species))
924 # take 0K energy where available, else total energy
925 if energy_0K:
926 energy = energy_0K
927 else:
928 energy = energy_total
929 # generate a minimal single-point calculator
930 sp_calc = SinglePointCalculator(atoms=atoms,
931 energy=energy,
932 forces=None,
933 magmoms=None,
934 stress=None)
935 atoms.calc = sp_calc
936 traj.append(atoms)
937 if index is None:
938 return traj
939 else:
940 return traj[index]
943def read_geom(filename, index=':', units=units_CODATA2002):
944 """
945 Wrapper function for the more generic read() functionality.
947 Note that this is function is intended to maintain backwards-compatibility
948 only. Keyword arguments will be passed to read_castep_geom().
949 """
950 from ase.io import read
951 return read(filename, index=index, format='castep-geom', units=units)
954def read_castep_geom(fd, index=None, units=units_CODATA2002):
955 """Reads a .geom file produced by the CASTEP GeometryOptimization task and
956 returns an atoms object.
957 The information about total free energy and forces of each atom for every
958 relaxation step will be stored for further analysis especially in a
959 single-point calculator.
960 Note that everything in the .geom file is in atomic units, which has
961 been conversed to commonly used unit angstrom(length) and eV (energy).
963 Note that the index argument has no effect as of now.
965 Contribution by Wei-Bing Zhang. Thanks!
967 Routine now accepts a filedescriptor in order to out-source the gz and
968 bz2 handling to formats.py. Note that there is a fallback routine
969 read_geom() that behaves like previous versions did.
970 """
971 from ase.calculators.singlepoint import SinglePointCalculator
973 # fd is closed by embracing read() routine
974 txt = fd.readlines()
976 traj = []
978 Hartree = units['Eh']
979 Bohr = units['a0']
981 # Yeah, we know that...
982 # print('N.B.: Energy in .geom file is not 0K extrapolated.')
983 for i, line in enumerate(txt):
984 if line.find('<-- E') > 0:
985 start_found = True
986 energy = float(line.split()[0]) * Hartree
987 cell = [x.split()[0:3] for x in txt[i + 1:i + 4]]
988 cell = np.array([[float(col) * Bohr for col in row] for row in
989 cell])
990 if line.find('<-- R') > 0 and start_found:
991 start_found = False
992 geom_start = i
993 for i, line in enumerate(txt[geom_start:]):
994 if line.find('<-- F') > 0:
995 geom_stop = i + geom_start
996 break
997 species = [line.split()[0] for line in
998 txt[geom_start:geom_stop]]
999 geom = np.array([[float(col) * Bohr for col in
1000 line.split()[2:5]] for line in
1001 txt[geom_start:geom_stop]])
1002 forces = np.array([[float(col) * Hartree / Bohr for col in
1003 line.split()[2:5]] for line in
1004 txt[geom_stop:geom_stop
1005 + (geom_stop - geom_start)]])
1006 image = ase.Atoms(species, geom, cell=cell, pbc=True)
1007 image.calc = SinglePointCalculator(
1008 atoms=image, energy=energy, forces=forces)
1009 traj.append(image)
1011 if index is None:
1012 return traj
1013 else:
1014 return traj[index]
1017def read_phonon(filename, index=None, read_vib_data=False,
1018 gamma_only=True, frequency_factor=None,
1019 units=units_CODATA2002):
1020 """
1021 Wrapper function for the more generic read() functionality.
1023 Note that this is function is intended to maintain backwards-compatibility
1024 only. For documentation see read_castep_phonon().
1025 """
1026 from ase.io import read
1028 if read_vib_data:
1029 full_output = True
1030 else:
1031 full_output = False
1033 return read(filename, index=index, format='castep-phonon',
1034 full_output=full_output, read_vib_data=read_vib_data,
1035 gamma_only=gamma_only, frequency_factor=frequency_factor,
1036 units=units)
1039def read_castep_phonon(fd, index=None, read_vib_data=False,
1040 gamma_only=True, frequency_factor=None,
1041 units=units_CODATA2002):
1042 """
1043 Reads a .phonon file written by a CASTEP Phonon task and returns an atoms
1044 object, as well as the calculated vibrational data if requested.
1046 Note that the index argument has no effect as of now.
1047 """
1049 # fd is closed by embracing read() routine
1050 lines = fd.readlines()
1052 atoms = None
1053 cell = []
1054 N = Nb = Nq = 0
1055 scaled_positions = []
1056 symbols = []
1057 masses = []
1059 # header
1060 L = 0
1061 while L < len(lines):
1063 line = lines[L]
1065 if 'Number of ions' in line:
1066 N = int(line.split()[3])
1067 elif 'Number of branches' in line:
1068 Nb = int(line.split()[3])
1069 elif 'Number of wavevectors' in line:
1070 Nq = int(line.split()[3])
1071 elif 'Unit cell vectors (A)' in line:
1072 for ll in range(3):
1073 L += 1
1074 fields = lines[L].split()
1075 cell.append([float(x) for x in fields[0:3]])
1076 elif 'Fractional Co-ordinates' in line:
1077 for ll in range(N):
1078 L += 1
1079 fields = lines[L].split()
1080 scaled_positions.append([float(x) for x in fields[1:4]])
1081 symbols.append(fields[4])
1082 masses.append(float(fields[5]))
1083 elif 'END header' in line:
1084 L += 1
1085 atoms = ase.Atoms(symbols=symbols,
1086 scaled_positions=scaled_positions,
1087 cell=cell)
1088 break
1090 L += 1
1092 # Eigenmodes and -vectors
1093 if frequency_factor is None:
1094 Kayser_to_eV = 1E2 * 2 * np.pi * units['hbar'] * units['c']
1095 # N.B. "fixed default" unit for frequencies in .phonon files is "cm-1"
1096 # (i.e. the latter is unaffected by the internal unit conversion system of
1097 # CASTEP!) set conversion factor to convert therefrom to eV by default for
1098 # now
1099 frequency_factor = Kayser_to_eV
1100 qpoints = []
1101 weights = []
1102 frequencies = []
1103 displacements = []
1104 for nq in range(Nq):
1105 fields = lines[L].split()
1106 qpoints.append([float(x) for x in fields[2:5]])
1107 weights.append(float(fields[5]))
1108 freqs = []
1109 for ll in range(Nb):
1110 L += 1
1111 fields = lines[L].split()
1112 freqs.append(frequency_factor * float(fields[1]))
1113 frequencies.append(np.array(freqs))
1115 # skip the two Phonon Eigenvectors header lines
1116 L += 2
1118 # generate a list of displacements with a structure that is identical to
1119 # what is stored internally in the Vibrations class (see in
1120 # ase.vibrations.Vibrations.modes):
1121 # np.array(displacements).shape == (Nb,3*N)
1123 disps = []
1124 for ll in range(Nb):
1125 disp_coords = []
1126 for lll in range(N):
1127 L += 1
1128 fields = lines[L].split()
1129 disp_x = float(fields[2]) + float(fields[3]) * 1.0j
1130 disp_y = float(fields[4]) + float(fields[5]) * 1.0j
1131 disp_z = float(fields[6]) + float(fields[7]) * 1.0j
1132 disp_coords.extend([disp_x, disp_y, disp_z])
1133 disps.append(np.array(disp_coords))
1134 displacements.append(np.array(disps))
1136 if read_vib_data:
1137 if gamma_only:
1138 vibdata = [frequencies[0], displacements[0]]
1139 else:
1140 vibdata = [qpoints, weights, frequencies, displacements]
1141 return vibdata, atoms
1142 else:
1143 return atoms
1146def read_md(filename, index=None, return_scalars=False,
1147 units=units_CODATA2002):
1148 """Wrapper function for the more generic read() functionality.
1150 Note that this function is intended to maintain backwards-compatibility
1151 only. For documentation see read_castep_md()
1152 """
1153 if return_scalars:
1154 full_output = True
1155 else:
1156 full_output = False
1158 from ase.io import read
1159 return read(filename, index=index, format='castep-md',
1160 full_output=full_output, return_scalars=return_scalars,
1161 units=units)
1164def read_castep_md(fd, index=None, return_scalars=False,
1165 units=units_CODATA2002):
1166 """Reads a .md file written by a CASTEP MolecularDynamics task
1167 and returns the trajectory stored therein as a list of atoms object.
1169 Note that the index argument has no effect as of now."""
1171 from ase.calculators.singlepoint import SinglePointCalculator
1173 factors = {
1174 't': units['t0'] * 1E15, # fs
1175 'E': units['Eh'], # eV
1176 'T': units['Eh'] / units['kB'],
1177 'P': units['Eh'] / units['a0']**3 * units['Pascal'],
1178 'h': units['a0'],
1179 'hv': units['a0'] / units['t0'],
1180 'S': units['Eh'] / units['a0']**3,
1181 'R': units['a0'],
1182 'V': np.sqrt(units['Eh'] / units['me']),
1183 'F': units['Eh'] / units['a0']}
1185 # fd is closed by embracing read() routine
1186 lines = fd.readlines()
1188 L = 0
1189 while 'END header' not in lines[L]:
1190 L += 1
1191 l_end_header = L
1192 lines = lines[l_end_header + 1:]
1193 times = []
1194 energies = []
1195 temperatures = []
1196 pressures = []
1197 traj = []
1199 # Initialization
1200 time = None
1201 Epot = None
1202 Ekin = None
1203 EH = None
1204 temperature = None
1205 pressure = None
1206 symbols = None
1207 positions = None
1208 cell = None
1209 velocities = None
1210 symbols = []
1211 positions = []
1212 velocities = []
1213 forces = []
1214 cell = np.eye(3)
1215 cell_velocities = []
1216 stress = []
1218 for (L, line) in enumerate(lines):
1219 fields = line.split()
1220 if len(fields) == 0:
1221 if L != 0:
1222 times.append(time)
1223 energies.append([Epot, EH, Ekin])
1224 temperatures.append(temperature)
1225 pressures.append(pressure)
1226 atoms = ase.Atoms(symbols=symbols,
1227 positions=positions,
1228 cell=cell)
1229 atoms.set_velocities(velocities)
1230 if len(stress) == 0:
1231 atoms.calc = SinglePointCalculator(
1232 atoms=atoms, energy=Epot, forces=forces)
1233 else:
1234 atoms.calc = SinglePointCalculator(
1235 atoms=atoms, energy=Epot,
1236 forces=forces, stress=stress)
1237 traj.append(atoms)
1238 symbols = []
1239 positions = []
1240 velocities = []
1241 forces = []
1242 cell = []
1243 cell_velocities = []
1244 stress = []
1245 continue
1246 if len(fields) == 1:
1247 time = factors['t'] * float(fields[0])
1248 continue
1250 if fields[-1] == 'E':
1251 E = [float(x) for x in fields[0:3]]
1252 Epot, EH, Ekin = [factors['E'] * Ei for Ei in E]
1253 continue
1255 if fields[-1] == 'T':
1256 temperature = factors['T'] * float(fields[0])
1257 continue
1259 # only printed in case of variable cell calculation or calculate_stress
1260 # explicitly requested
1261 if fields[-1] == 'P':
1262 pressure = factors['P'] * float(fields[0])
1263 continue
1264 if fields[-1] == 'h':
1265 h = [float(x) for x in fields[0:3]]
1266 cell.append([factors['h'] * hi for hi in h])
1267 continue
1269 # only printed in case of variable cell calculation
1270 if fields[-1] == 'hv':
1271 hv = [float(x) for x in fields[0:3]]
1272 cell_velocities.append([factors['hv'] * hvi for hvi in hv])
1273 continue
1275 # only printed in case of variable cell calculation
1276 if fields[-1] == 'S':
1277 S = [float(x) for x in fields[0:3]]
1278 stress.append([factors['S'] * Si for Si in S])
1279 continue
1280 if fields[-1] == 'R':
1281 symbols.append(fields[0])
1282 R = [float(x) for x in fields[2:5]]
1283 positions.append([factors['R'] * Ri for Ri in R])
1284 continue
1285 if fields[-1] == 'V':
1286 V = [float(x) for x in fields[2:5]]
1287 velocities.append([factors['V'] * Vi for Vi in V])
1288 continue
1289 if fields[-1] == 'F':
1290 F = [float(x) for x in fields[2:5]]
1291 forces.append([factors['F'] * Fi for Fi in F])
1292 continue
1294 if index is None:
1295 pass
1296 else:
1297 traj = traj[index]
1299 if return_scalars:
1300 data = [times, energies, temperatures, pressures]
1301 return data, traj
1302 else:
1303 return traj
1306# Routines that only the calculator requires
1308def read_param(filename='', calc=None, fd=None, get_interface_options=False):
1309 if fd is None:
1310 if filename == '':
1311 raise ValueError('One between filename and fd must be provided')
1312 fd = open(filename)
1313 elif filename:
1314 warnings.warn('Filestream used to read param, file name will be '
1315 'ignored')
1317 # If necessary, get the interface options
1318 if get_interface_options:
1319 int_opts = {}
1320 optre = re.compile(r'# ASE_INTERFACE ([^\s]+) : ([^\s]+)')
1322 lines = fd.readlines()
1323 fd.seek(0)
1325 for line in lines:
1326 m = optre.search(line)
1327 if m:
1328 int_opts[m.groups()[0]] = m.groups()[1]
1330 data = read_freeform(fd)
1332 if calc is None:
1333 from ase.calculators.castep import Castep
1334 calc = Castep(check_castep_version=False, keyword_tolerance=2)
1336 for kw, (val, otype) in data.items():
1337 if otype == 'block':
1338 val = val.split('\n') # Avoids a bug for one-line blocks
1339 calc.param.__setattr__(kw, val)
1341 if not get_interface_options:
1342 return calc
1343 else:
1344 return calc, int_opts
1347def write_param(filename, param, check_checkfile=False,
1348 force_write=False,
1349 interface_options=None):
1350 """Writes a CastepParam object to a CASTEP .param file
1352 Parameters:
1353 filename: the location of the file to write to. If it
1354 exists it will be overwritten without warning. If it
1355 doesn't it will be created.
1356 param: a CastepParam instance
1357 check_checkfile : if set to True, write_param will
1358 only write continuation or reuse statement
1359 if a restart file exists in the same directory
1360 """
1361 if os.path.isfile(filename) and not force_write:
1362 warnings.warn('ase.io.castep.write_param: Set optional argument '
1363 'force_write=True to overwrite %s.' % filename)
1364 return False
1366 out = paropen(filename, 'w')
1367 out.write('#######################################################\n')
1368 out.write('#CASTEP param file: %s\n' % filename)
1369 out.write('#Created using the Atomic Simulation Environment (ASE)#\n')
1370 if interface_options is not None:
1371 out.write('# Internal settings of the calculator\n')
1372 out.write('# This can be switched off by settings\n')
1373 out.write('# calc._export_settings = False\n')
1374 out.write('# If stated, this will be automatically processed\n')
1375 out.write('# by ase.io.castep.read_seed()\n')
1376 for option, value in sorted(interface_options.items()):
1377 out.write('# ASE_INTERFACE %s : %s\n' % (option, value))
1378 out.write('#######################################################\n\n')
1380 if check_checkfile:
1381 param = deepcopy(param) # To avoid modifying the parent one
1382 for checktype in ['continuation', 'reuse']:
1383 opt = getattr(param, checktype)
1384 if opt and opt.value:
1385 fname = opt.value
1386 if fname == 'default':
1387 fname = os.path.splitext(filename)[0] + '.check'
1388 if not (os.path.exists(fname) or
1389 # CASTEP also understands relative path names, hence
1390 # also check relative to the param file directory
1391 os.path.exists(
1392 os.path.join(os.path.dirname(filename),
1393 opt.value))):
1394 opt.clear()
1396 write_freeform(out, param)
1398 out.close()
1401def read_seed(seed, new_seed=None, ignore_internal_keys=False):
1402 """A wrapper around the CASTEP Calculator in conjunction with
1403 read_cell and read_param. Basically this can be used to reuse
1404 a previous calculation which results in a triple of
1405 cell/param/castep file. The label of the calculation if pre-
1406 fixed with `copy_of_` and everything else will be recycled as
1407 much as possible from the addressed calculation.
1409 Please note that this routine will return an atoms ordering as specified
1410 in the cell file! It will thus undo the potential reordering internally
1411 done by castep.
1412 """
1414 directory = os.path.abspath(os.path.dirname(seed))
1415 seed = os.path.basename(seed)
1417 paramfile = os.path.join(directory, '%s.param' % seed)
1418 cellfile = os.path.join(directory, '%s.cell' % seed)
1419 castepfile = os.path.join(directory, '%s.castep' % seed)
1420 checkfile = os.path.join(directory, '%s.check' % seed)
1422 atoms = read_cell(cellfile)
1423 atoms.calc._directory = directory
1424 atoms.calc._rename_existing_dir = False
1425 atoms.calc._castep_pp_path = directory
1426 atoms.calc.merge_param(paramfile,
1427 ignore_internal_keys=ignore_internal_keys)
1428 if new_seed is None:
1429 atoms.calc._label = 'copy_of_%s' % seed
1430 else:
1431 atoms.calc._label = str(new_seed)
1432 if os.path.isfile(castepfile):
1433 # _set_atoms needs to be True here
1434 # but we set it right back to False
1435 # atoms.calc._set_atoms = False
1436 # BUGFIX: I do not see a reason to do that!
1437 atoms.calc.read(castepfile)
1438 # atoms.calc._set_atoms = False
1440 # if here is a check file, we also want to re-use this information
1441 if os.path.isfile(checkfile):
1442 atoms.calc._check_file = os.path.basename(checkfile)
1444 # sync the top-level object with the
1445 # one attached to the calculator
1446 atoms = atoms.calc.atoms
1447 else:
1448 # There are cases where we only want to restore a calculator/atoms
1449 # setting without a castep file...
1450 pass
1451 # No print statement required in these cases
1452 warnings.warn(
1453 'Corresponding *.castep file not found. '
1454 'Atoms object will be restored from *.cell and *.param only.')
1455 atoms.calc.push_oldstate()
1457 return atoms
1460def read_bands(filename='', fd=None, units=units_CODATA2002):
1461 """Read Castep.bands file to kpoints, weights and eigenvalues
1463 Args:
1464 filename (str):
1465 path to seedname.bands file
1466 fd (fd):
1467 file descriptor for open bands file
1468 units (dict):
1469 Conversion factors for atomic units
1471 Returns:
1472 (tuple):
1473 (kpts, weights, eigenvalues, efermi)
1475 Where ``kpts`` and ``weights`` are 1d numpy arrays, eigenvalues
1476 is an array of shape (spin, kpts, nbands) and efermi is a float
1477 """
1479 Hartree = units['Eh']
1481 if fd is None:
1482 if filename == '':
1483 raise ValueError('One between filename and fd must be provided')
1484 fd = open(filename)
1485 elif filename:
1486 warnings.warn('Filestream used to read param, file name will be '
1487 'ignored')
1489 nkpts, nspin, _, nbands, efermi = [t(fd.readline().split()[-1]) for t in
1490 [int, int, float, int, float]]
1492 kpts, weights = np.zeros((nkpts, 3)), np.zeros(nkpts)
1493 eigenvalues = np.zeros((nspin, nkpts, nbands))
1495 # Skip unit cell
1496 for _ in range(4):
1497 fd.readline()
1499 def _kptline_to_i_k_wt(line):
1500 line = line.split()
1501 line = [int(line[1])] + list(map(float, line[2:]))
1502 return (line[0] - 1, line[1:4], line[4])
1504 # CASTEP often writes these out-of-order, so check index and write directly
1505 # to the correct row
1506 for kpt_line in range(nkpts):
1507 i_kpt, kpt, wt = _kptline_to_i_k_wt(fd.readline())
1508 kpts[i_kpt, :], weights[i_kpt] = kpt, wt
1509 for spin in range(nspin):
1510 fd.readline() # Skip 'Spin component N' line
1511 eigenvalues[spin, i_kpt, :] = [float(fd.readline())
1512 for _ in range(nbands)]
1514 return (kpts, weights, eigenvalues * Hartree, efermi * Hartree)