Coverage for ase / io / castep / __init__.py: 74.23%
555 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-21 15:52 +0000
1# fmt: off
3"""This module defines I/O routines with CASTEP files.
4The key idea is that all function accept or return atoms objects.
5CASTEP specific parameters will be returned through the <atoms>.calc
6attribute.
7"""
8import os
9import re
10import warnings
11from copy import deepcopy
13import numpy as np
15import ase
17# independent unit management included here:
18# When high accuracy is required, this allows to easily pin down
19# unit conversion factors from different "unit definition systems"
20# (CODATA1986 for ase-3.6.0.2515 vs CODATA2002 for CASTEP 5.01).
21#
22# ase.units in in ase-3.6.0.2515 is based on CODATA1986
23import ase.units
24from ase.constraints import FixAtoms, FixCartesian, FixedLine, FixedPlane
25from ase.geometry.cell import cellpar_to_cell
26from ase.io.castep.castep_reader import read_castep_castep
27from ase.parallel import paropen
28from ase.spacegroup import Spacegroup
29from ase.utils import (
30 atoms_to_spglib_cell,
31 reader,
32 spglib_new_errorhandling,
33 writer,
34)
36from .geom_md_ts import (
37 read_castep_geom,
38 read_castep_md,
39 write_castep_geom,
40 write_castep_md,
41)
43units_ase = {
44 'hbar': ase.units._hbar * ase.units.J,
45 'Eh': ase.units.Hartree,
46 'kB': ase.units.kB,
47 'a0': ase.units.Bohr,
48 't0': ase.units._hbar * ase.units.J / ase.units.Hartree,
49 'c': ase.units._c,
50 'me': ase.units._me / ase.units._amu,
51 'Pascal': 1.0 / ase.units.Pascal}
53# CODATA1986 (included herein for the sake of completeness)
54# taken from
55# http://physics.nist.gov/cuu/Archive/1986RMP.pdf
56units_CODATA1986 = {
57 'hbar': 6.5821220E-16, # eVs
58 'Eh': 27.2113961, # eV
59 'kB': 8.617385E-5, # eV/K
60 'a0': 0.529177249, # A
61 'c': 299792458, # m/s
62 'e': 1.60217733E-19, # C
63 'me': 5.485799110E-4} # u
65# CODATA2002: default in CASTEP 5.01
66# (-> check in more recent CASTEP in case of numerical discrepancies?!)
67# taken from
68# http://physics.nist.gov/cuu/Document/all_2002.pdf
69units_CODATA2002 = {
70 'hbar': 6.58211915E-16, # eVs
71 'Eh': 27.2113845, # eV
72 'kB': 8.617343E-5, # eV/K
73 'a0': 0.5291772108, # A
74 'c': 299792458, # m/s
75 'e': 1.60217653E-19, # C
76 'me': 5.4857990945E-4} # u
78# (common) derived entries
79for d in (units_CODATA1986, units_CODATA2002):
80 d['t0'] = d['hbar'] / d['Eh'] # s
81 d['Pascal'] = d['e'] * 1E30 # Pa
84__all__ = [
85 # routines for the generic io function
86 'read_castep_castep',
87 'read_castep_cell',
88 'read_castep_geom',
89 'read_castep_md',
90 'read_phonon',
91 'read_castep_phonon',
92 # additional reads that still need to be wrapped
93 'read_param',
94 'read_seed',
95 # write that is already wrapped
96 'write_castep_cell',
97 'write_castep_geom',
98 'write_castep_md',
99 # param write - in principle only necessary in junction with the calculator
100 'write_param']
103def write_freeform(fd, outputobj):
104 """
105 Prints out to a given file a CastepInputFile or derived class, such as
106 CastepCell or CastepParam.
107 """
109 options = outputobj._options
111 # Some keywords, if present, are printed in this order
112 preferred_order = ['lattice_cart', 'lattice_abc',
113 'positions_frac', 'positions_abs',
114 'species_pot', 'symmetry_ops', # CELL file
115 'task', 'cut_off_energy' # PARAM file
116 ]
118 keys = outputobj.get_attr_dict().keys()
119 # This sorts only the ones in preferred_order and leaves the rest
120 # untouched
121 keys = sorted(keys, key=lambda x: preferred_order.index(x)
122 if x in preferred_order
123 else len(preferred_order))
125 for kw in keys:
126 opt = options[kw]
127 if opt.type.lower() == 'block':
128 fd.write('%BLOCK {0}\n{1}\n%ENDBLOCK {0}\n\n'.format(
129 kw.upper(),
130 opt.value.strip('\n')))
131 else:
132 fd.write(f'{kw.upper()}: {opt.value}\n')
135@writer
136def write_castep_cell(fd, atoms, positions_frac=False, force_write=False,
137 precision=6, magnetic_moments=None,
138 castep_cell=None):
139 """
140 This CASTEP export function write minimal information to
141 a .cell file. If the atoms object is a trajectory, it will
142 take the last image.
144 Note that function has been altered in order to require a filedescriptor
145 rather than a filename. This allows to use the more generic write()
146 function from formats.py
148 Note that the "force_write" keywords has no effect currently.
150 Arguments:
152 positions_frac: boolean. If true, positions are printed as fractional
153 rather than absolute. Default is false.
154 castep_cell: if provided, overrides the existing CastepCell object in
155 the Atoms calculator
156 precision: number of digits to which lattice and positions are printed
157 magnetic_moments: if None, no SPIN values are initialised.
158 If 'initial', the values from
159 get_initial_magnetic_moments() are used.
160 If 'calculated', the values from
161 get_magnetic_moments() are used.
162 If an array of the same length as the atoms object,
163 its contents will be used as magnetic moments.
164 """
166 if isinstance(atoms, list):
167 if len(atoms) > 1:
168 atoms = atoms[-1]
170 # Header
171 fd.write('# written by ASE\n\n')
173 # To write this we simply use the existing Castep calculator, or create
174 # one
175 from ase.calculators.castep import Castep, CastepCell
177 try:
178 has_cell = isinstance(atoms.calc.cell, CastepCell)
179 except AttributeError:
180 has_cell = False
182 if has_cell:
183 cell = deepcopy(atoms.calc.cell)
184 else:
185 cell = Castep(keyword_tolerance=2).cell
187 # Write lattice
188 fformat = f'%{precision + 3}.{precision}f'
189 cell_block_format = ' '.join([fformat] * 3)
190 cell.lattice_cart = [cell_block_format % tuple(line)
191 for line in atoms.get_cell()]
193 if positions_frac:
194 pos_keyword = 'positions_frac'
195 positions = atoms.get_scaled_positions()
196 else:
197 pos_keyword = 'positions_abs'
198 positions = atoms.get_positions()
200 if atoms.has('castep_custom_species'):
201 elems = atoms.get_array('castep_custom_species')
202 else:
203 elems = atoms.get_chemical_symbols()
204 if atoms.has('masses'):
206 from ase.data import atomic_masses
207 masses = atoms.get_array('masses')
208 custom_masses = {}
210 for i, species in enumerate(elems):
211 custom_mass = masses[i]
213 # build record of different masses for each species
214 if species not in custom_masses:
216 # build dictionary of positions of all species with
217 # same name and mass value ideally there should only
218 # be one mass per species
219 custom_masses[species] = {custom_mass: [i]}
221 # if multiple masses found for a species
222 elif custom_mass not in custom_masses[species].keys():
224 # if custom species were already manually defined raise an error
225 if atoms.has('castep_custom_species'):
226 raise ValueError(
227 "Could not write custom mass block for {0}. \n"
228 "Custom mass was set ({1}), but an inconsistent set of "
229 "castep_custom_species already defines "
230 "({2}) for {0}. \n"
231 "If using both features, ensure that "
232 "each species type in "
233 "atoms.arrays['castep_custom_species'] "
234 "has consistent mass values and that each atom "
235 "with non-standard "
236 "mass belongs to a custom species type."
237 "".format(
238 species, custom_mass, list(
239 custom_masses[species].keys())[0]))
241 # append mass to create custom species later
242 else:
243 custom_masses[species][custom_mass] = [i]
244 else:
245 custom_masses[species][custom_mass].append(i)
247 # create species_mass block
248 mass_block = []
250 for el, mass_dict in custom_masses.items():
252 # ignore mass record that match defaults
253 default = mass_dict.pop(atomic_masses[atoms.get_array(
254 'numbers')[list(elems).index(el)]], None)
255 if mass_dict:
256 # no custom species need to be created
257 if len(mass_dict) == 1 and not default:
258 mass_block.append('{} {}'.format(
259 el, list(mass_dict.keys())[0]))
260 # for each custom mass, create new species and change names to
261 # match in 'elems' list
262 else:
263 warnings.warn(
264 'Custom mass specified for '
265 'standard species {}, creating custom species'
266 .format(el))
268 for i, vals in enumerate(mass_dict.items()):
269 mass_val, idxs = vals
270 custom_species_name = f"{el}:{i}"
271 warnings.warn(
272 'Creating custom species {} with mass {}'.format(
273 custom_species_name, str(mass_dict)))
274 for idx in idxs:
275 elems[idx] = custom_species_name
276 mass_block.append('{} {}'.format(
277 custom_species_name, mass_val))
279 cell.species_mass = mass_block
281 if atoms.has('castep_labels'):
282 labels = atoms.get_array('castep_labels')
283 else:
284 labels = ['NULL'] * len(elems)
286 if str(magnetic_moments).lower() == 'initial':
287 magmoms = atoms.get_initial_magnetic_moments()
288 elif str(magnetic_moments).lower() == 'calculated':
289 magmoms = atoms.get_magnetic_moments()
290 elif np.array(magnetic_moments).shape == (len(elems),):
291 magmoms = np.array(magnetic_moments)
292 else:
293 magmoms = [0] * len(elems)
295 pos_block = []
296 pos_block_format = '%s ' + cell_block_format
298 for i, el in enumerate(elems):
299 xyz = positions[i]
300 line = pos_block_format % tuple([el] + list(xyz))
301 # ADD other keywords if necessary
302 if magmoms[i] != 0:
303 line += f' SPIN={magmoms[i]} '
304 if labels[i].strip() not in ('NULL', ''):
305 line += f' LABEL={labels[i]} '
306 pos_block.append(line)
308 setattr(cell, pos_keyword, pos_block)
310 constr_block = _make_block_ionic_constraints(atoms)
311 if constr_block:
312 cell.ionic_constraints = constr_block
314 write_freeform(fd, cell)
317def _make_block_ionic_constraints(atoms: ase.Atoms) -> list[str]:
318 constr_block: list[str] = []
319 species_indices = atoms.symbols.species_indices()
320 for constr in atoms.constraints:
321 if not _is_constraint_valid(constr, len(atoms)):
322 continue
323 for i in constr.index:
324 symbol = atoms.get_chemical_symbols()[i]
325 nis = species_indices[i] + 1
326 if isinstance(constr, FixAtoms):
327 for j in range(3): # constraint for all three directions
328 ic = len(constr_block) + 1
329 line = f'{ic:6d} {symbol:3s} {nis:3d} '
330 line += ['1 0 0', '0 1 0', '0 0 1'][j]
331 constr_block.append(line)
332 elif isinstance(constr, FixCartesian):
333 for j, m in enumerate(constr.mask):
334 if m == 0: # not constrained
335 continue
336 ic = len(constr_block) + 1
337 line = f'{ic:6d} {symbol:3s} {nis:3d} '
338 line += ['1 0 0', '0 1 0', '0 0 1'][j]
339 constr_block.append(line)
340 elif isinstance(constr, FixedPlane):
341 ic = len(constr_block) + 1
342 line = f'{ic:6d} {symbol:3s} {nis:3d} '
343 line += ' '.join([str(d) for d in constr.dir])
344 constr_block.append(line)
345 elif isinstance(constr, FixedLine):
346 for direction in _calc_normal_vectors(constr):
347 ic = len(constr_block) + 1
348 line = f'{ic:6d} {symbol:3s} {nis:3d} '
349 line += ' '.join(str(_) for _ in direction)
350 constr_block.append(line)
351 return constr_block
354def _is_constraint_valid(constraint, natoms: int) -> bool:
355 supported_constraints = (FixAtoms, FixedPlane, FixedLine, FixCartesian)
356 if not isinstance(constraint, supported_constraints):
357 warnings.warn(f'{constraint} is not supported by ASE CASTEP, skipped')
358 return False
359 if any(i < 0 or i >= natoms for i in constraint.index):
360 warnings.warn(f'{constraint} contains invalid indices, skipped')
361 return False
362 return True
365def _calc_normal_vectors(constr: FixedLine) -> tuple[np.ndarray, np.ndarray]:
366 direction = constr.dir
368 i2, i1 = np.argsort(np.abs(direction))[1:]
369 v1 = direction[i1]
370 v2 = direction[i2]
371 n1 = np.zeros(3)
372 n1[i2] = v1
373 n1[i1] = -v2
374 n1 = n1 / np.linalg.norm(n1)
376 n2 = np.cross(direction, n1)
377 n2 = n2 / np.linalg.norm(n2)
379 return n1, n2
382def read_freeform(fd):
383 """
384 Read a CASTEP freeform file (the basic format of .cell and .param files)
385 and return keyword-value pairs as a dict (values are strings for single
386 keywords and lists of strings for blocks).
387 """
389 from ase.io.castep.castep_input_file import CastepInputFile
391 inputobj = CastepInputFile(keyword_tolerance=2)
393 filelines = fd.readlines()
395 keyw = None
396 read_block = False
397 block_lines = None
399 for i, l in enumerate(filelines):
401 # Strip all comments, aka anything after a hash
402 L = re.split(r'[#!;]', l, maxsplit=1)[0].strip()
404 if L == '':
405 # Empty line... skip
406 continue
408 lsplit = re.split(r'\s*[:=]*\s+', L, maxsplit=1)
410 if read_block:
411 if lsplit[0].lower() == '%endblock':
412 if len(lsplit) == 1 or lsplit[1].lower() != keyw:
413 raise ValueError('Out of place end of block at '
414 'line %i in freeform file' % i + 1)
415 else:
416 read_block = False
417 inputobj.__setattr__(keyw, block_lines)
418 else:
419 block_lines += [L]
420 else:
421 # Check the first word
423 # Is it a block?
424 read_block = (lsplit[0].lower() == '%block')
425 if read_block:
426 if len(lsplit) == 1:
427 raise ValueError(('Unrecognizable block at line %i '
428 'in io freeform file') % i + 1)
429 else:
430 keyw = lsplit[1].lower()
431 else:
432 keyw = lsplit[0].lower()
434 # Now save the value
435 if read_block:
436 block_lines = []
437 else:
438 inputobj.__setattr__(keyw, ' '.join(lsplit[1:]))
440 return inputobj.get_attr_dict(types=True)
443@reader
444def read_castep_cell(fd, index=None, calculator_args={}, find_spg=False,
445 units=units_CODATA2002):
446 """Read a .cell file and return an atoms object.
447 Any value found that does not fit the atoms API
448 will be stored in the atoms.calc attribute.
450 By default, the Castep calculator will be tolerant and in the absence of a
451 castep_keywords.json file it will just accept all keywords that aren't
452 automatically parsed.
453 """
455 from ase.calculators.castep import Castep
457 cell_units = { # Units specifiers for CASTEP
458 'bohr': units_CODATA2002['a0'],
459 'ang': 1.0,
460 'm': 1e10,
461 'cm': 1e8,
462 'nm': 10,
463 'pm': 1e-2
464 }
466 calc = Castep(**calculator_args)
468 if calc.cell.castep_version == 0 and calc._kw_tol < 3:
469 # No valid castep_keywords.json was found
470 warnings.warn(
471 'read_cell: Warning - Was not able to validate CASTEP input. '
472 'This may be due to a non-existing '
473 '"castep_keywords.json" '
474 'file or a non-existing CASTEP installation. '
475 'Parsing will go on but keywords will not be '
476 'validated and may cause problems if incorrect during a CASTEP '
477 'run.')
479 celldict = read_freeform(fd)
481 def parse_blockunit(line_tokens, blockname):
482 u = 1.0
483 if len(line_tokens[0]) == 1:
484 usymb = line_tokens[0][0].lower()
485 u = cell_units.get(usymb, 1)
486 if usymb not in cell_units:
487 warnings.warn('read_cell: Warning - ignoring invalid '
488 'unit specifier in %BLOCK {} '
489 '(assuming Angstrom instead)'.format(blockname))
490 line_tokens = line_tokens[1:]
491 return u, line_tokens
493 # Arguments to pass to the Atoms object at the end
494 aargs = {
495 'pbc': True
496 }
498 # Start by looking for the lattice
499 lat_keywords = [w in celldict for w in ('lattice_cart', 'lattice_abc')]
500 if all(lat_keywords):
501 warnings.warn('read_cell: Warning - two lattice blocks present in the'
502 ' same file. LATTICE_ABC will be ignored')
503 elif not any(lat_keywords):
504 raise ValueError('Cell file must contain at least one between '
505 'LATTICE_ABC and LATTICE_CART')
507 if 'lattice_abc' in celldict:
509 lines = celldict.pop('lattice_abc')[0].split('\n')
510 line_tokens = [line.split() for line in lines]
512 u, line_tokens = parse_blockunit(line_tokens, 'lattice_abc')
514 if len(line_tokens) != 2:
515 warnings.warn('read_cell: Warning - ignoring additional '
516 'lines in invalid %BLOCK LATTICE_ABC')
518 abc = [float(p) * u for p in line_tokens[0][:3]]
519 angles = [float(phi) for phi in line_tokens[1][:3]]
521 aargs['cell'] = cellpar_to_cell(abc + angles)
523 if 'lattice_cart' in celldict:
525 lines = celldict.pop('lattice_cart')[0].split('\n')
526 line_tokens = [line.split() for line in lines]
528 u, line_tokens = parse_blockunit(line_tokens, 'lattice_cart')
530 if len(line_tokens) != 3:
531 warnings.warn('read_cell: Warning - ignoring more than '
532 'three lattice vectors in invalid %BLOCK '
533 'LATTICE_CART')
535 aargs['cell'] = [[float(x) * u for x in lt[:3]] for lt in line_tokens]
537 # Now move on to the positions
538 pos_keywords = [w in celldict
539 for w in ('positions_abs', 'positions_frac')]
541 if all(pos_keywords):
542 warnings.warn('read_cell: Warning - two lattice blocks present in the'
543 ' same file. POSITIONS_FRAC will be ignored')
544 del celldict['positions_frac']
545 elif not any(pos_keywords):
546 raise ValueError('Cell file must contain at least one between '
547 'POSITIONS_FRAC and POSITIONS_ABS')
549 aargs['symbols'] = []
550 pos_type = 'positions'
551 pos_block = celldict.pop('positions_abs', [None])[0]
552 if pos_block is None:
553 pos_type = 'scaled_positions'
554 pos_block = celldict.pop('positions_frac', [None])[0]
555 aargs[pos_type] = []
557 lines = pos_block.split('\n')
558 line_tokens = [line.split() for line in lines]
560 if 'scaled' not in pos_type:
561 u, line_tokens = parse_blockunit(line_tokens, 'positions_abs')
562 else:
563 u = 1.0
565 # Here we extract all the possible additional info
566 # These are marked by their type
568 add_info = {
569 'SPIN': (float, 0.0), # (type, default)
570 'MAGMOM': (float, 0.0),
571 'LABEL': (str, 'NULL')
572 }
573 add_info_arrays = {k: [] for k in add_info}
575 def parse_info(raw_info):
577 re_keys = (r'({0})\s*[=:\s]{{1}}\s'
578 r'*([^\s]*)').format('|'.join(add_info.keys()))
579 # Capture all info groups
580 info = re.findall(re_keys, raw_info)
581 info = {g[0]: add_info[g[0]][0](g[1]) for g in info}
582 return info
584 # Array for custom species (a CASTEP special thing)
585 # Usually left unused
586 custom_species = None
588 for tokens in line_tokens:
589 # Now, process the whole 'species' thing
590 spec_custom = tokens[0].split(':', 1)
591 elem = spec_custom[0]
592 if len(spec_custom) > 1 and custom_species is None:
593 # Add it to the custom info!
594 custom_species = list(aargs['symbols'])
595 if custom_species is not None:
596 custom_species.append(tokens[0])
597 aargs['symbols'].append(elem)
598 aargs[pos_type].append([float(p) * u for p in tokens[1:4]])
599 # Now for the additional information
600 info = ' '.join(tokens[4:])
601 info = parse_info(info)
602 for k in add_info:
603 add_info_arrays[k] += [info.get(k, add_info[k][1])]
605 # read in custom species mass
606 if 'species_mass' in celldict:
607 spec_list = custom_species if custom_species else aargs['symbols']
608 aargs['masses'] = [None for _ in spec_list]
609 lines = celldict.pop('species_mass')[0].split('\n')
610 line_tokens = [line.split() for line in lines]
612 if len(line_tokens[0]) == 1:
613 if line_tokens[0][0].lower() not in ('amu', 'u'):
614 raise ValueError(
615 "unit specifier '{}' in %BLOCK SPECIES_MASS "
616 "not recognised".format(
617 line_tokens[0][0].lower()))
618 line_tokens = line_tokens[1:]
620 for tokens in line_tokens:
621 token_pos_list = [i for i, x in enumerate(
622 spec_list) if x == tokens[0]]
623 if len(token_pos_list) == 0:
624 warnings.warn(
625 'read_cell: Warning - ignoring unused '
626 'species mass {} in %BLOCK SPECIES_MASS'.format(
627 tokens[0]))
628 for idx in token_pos_list:
629 aargs['masses'][idx] = tokens[1]
631 # Now on to the species potentials...
632 if 'species_pot' in celldict:
633 lines = celldict.pop('species_pot')[0].split('\n')
634 line_tokens = [line.split() for line in lines]
636 for tokens in line_tokens:
637 if len(tokens) == 1:
638 # It's a library
639 all_spec = (set(custom_species) if custom_species is not None
640 else set(aargs['symbols']))
641 for s in all_spec:
642 calc.cell.species_pot = (s, tokens[0])
643 else:
644 calc.cell.species_pot = tuple(tokens[:2])
646 # Ionic constraints
647 raw_constraints = {}
649 if 'ionic_constraints' in celldict:
650 lines = celldict.pop('ionic_constraints')[0].split('\n')
651 line_tokens = [line.split() for line in lines]
653 for tokens in line_tokens:
654 if len(tokens) != 6:
655 continue
656 _, species, nic, *xyz = tokens
658 nic = int(nic)
659 raw_constraints.setdefault((species, nic), [])
660 raw_constraints[(species, nic)].append(np.array(xyz, dtype=float))
662 # Symmetry operations
663 if 'symmetry_ops' in celldict:
664 lines = celldict.pop('symmetry_ops')[0].split('\n')
665 line_tokens = [line.split() for line in lines]
667 # Read them in blocks of four
668 blocks = np.array(line_tokens).astype(float)
669 if (len(blocks.shape) != 2 or blocks.shape[1] != 3
670 or blocks.shape[0] % 4 != 0):
671 warnings.warn('Warning: could not parse SYMMETRY_OPS'
672 ' block properly, skipping')
673 else:
674 blocks = blocks.reshape((-1, 4, 3))
675 rotations = blocks[:, :3]
676 translations = blocks[:, 3]
678 # Regardless of whether we recognize them, store these
679 calc.cell.symmetry_ops = (rotations, translations)
681 # Anything else that remains, just add it to the cell object:
682 for k, (val, otype) in celldict.items():
683 try:
684 if otype == 'block':
685 val = val.split('\n') # Avoids a bug for one-line blocks
686 calc.cell.__setattr__(k, val)
687 except Exception as e:
688 raise RuntimeError(
689 f'Problem setting calc.cell.{k} = {val}: {e}')
691 # Get the relevant additional info
692 aargs['magmoms'] = np.array(add_info_arrays['SPIN'])
693 # SPIN or MAGMOM are alternative keywords
694 aargs['magmoms'] = np.where(aargs['magmoms'] != 0,
695 aargs['magmoms'],
696 add_info_arrays['MAGMOM'])
697 labels = np.array(add_info_arrays['LABEL'])
699 aargs['calculator'] = calc
701 atoms = ase.Atoms(**aargs)
702 _set_spacegroup_info(find_spg, atoms)
704 atoms.new_array('castep_labels', labels)
705 if custom_species is not None:
706 atoms.new_array('castep_custom_species', np.array(custom_species))
708 fixed_atoms = []
709 constraints = []
710 index_dict = atoms.symbols.indices()
711 for (species, nic), value in raw_constraints.items():
713 absolute_nr = index_dict[species][nic - 1]
714 if len(value) == 3:
715 # Check if they are linearly independent
716 if np.linalg.det(value) == 0:
717 warnings.warn(
718 'Error: Found linearly dependent constraints attached '
719 'to atoms %s' %
720 (absolute_nr))
721 continue
722 fixed_atoms.append(absolute_nr)
723 elif len(value) == 2:
724 direction = np.cross(value[0], value[1])
725 # Check if they are linearly independent
726 if np.linalg.norm(direction) == 0:
727 warnings.warn(
728 'Error: Found linearly dependent constraints attached '
729 'to atoms %s' %
730 (absolute_nr))
731 continue
732 constraint = FixedLine(indices=absolute_nr, direction=direction)
733 constraints.append(constraint)
734 elif len(value) == 1:
735 direction = np.array(value[0], dtype=float)
736 constraint = FixedPlane(indices=absolute_nr, direction=direction)
737 constraints.append(constraint)
738 else:
739 warnings.warn(
740 f'Error: Found {len(value)} statements attached to atoms '
741 f'{absolute_nr}'
742 )
744 # we need to sort the fixed atoms list in order not to raise an assertion
745 # error in FixAtoms
746 if fixed_atoms:
747 constraints.append(FixAtoms(indices=sorted(fixed_atoms)))
748 if constraints:
749 atoms.set_constraint(constraints)
751 atoms.calc.atoms = atoms
752 atoms.calc.push_oldstate()
754 return atoms
757def read_phonon(filename, index=None, read_vib_data=False,
758 gamma_only=True, frequency_factor=None,
759 units=units_CODATA2002):
760 """
761 Wrapper function for the more generic read() functionality.
763 Note that this is function is intended to maintain backwards-compatibility
764 only. For documentation see read_castep_phonon().
765 """
766 from ase.io import read
768 if read_vib_data:
769 full_output = True
770 else:
771 full_output = False
773 return read(filename, index=index, format='castep-phonon',
774 full_output=full_output, read_vib_data=read_vib_data,
775 gamma_only=gamma_only, frequency_factor=frequency_factor,
776 units=units)
779def read_castep_phonon(fd, index=None, read_vib_data=False,
780 gamma_only=True, frequency_factor=None,
781 units=units_CODATA2002):
782 """
783 Reads a .phonon file written by a CASTEP Phonon task and returns an atoms
784 object, as well as the calculated vibrational data if requested.
786 Note that the index argument has no effect as of now.
787 """
789 # fd is closed by embracing read() routine
790 lines = fd.readlines()
792 atoms = None
793 cell = []
794 N = Nb = Nq = 0
795 scaled_positions = []
796 symbols = []
797 masses = []
799 # header
800 L = 0
801 while L < len(lines):
803 line = lines[L]
805 if 'Number of ions' in line:
806 N = int(line.split()[3])
807 elif 'Number of branches' in line:
808 Nb = int(line.split()[3])
809 elif 'Number of wavevectors' in line:
810 Nq = int(line.split()[3])
811 elif 'Unit cell vectors (A)' in line:
812 for _ in range(3):
813 L += 1
814 fields = lines[L].split()
815 cell.append([float(x) for x in fields[0:3]])
816 elif 'Fractional Co-ordinates' in line:
817 for _ in range(N):
818 L += 1
819 fields = lines[L].split()
820 scaled_positions.append([float(x) for x in fields[1:4]])
821 symbols.append(fields[4])
822 masses.append(float(fields[5]))
823 elif 'END header' in line:
824 L += 1
825 atoms = ase.Atoms(symbols=symbols,
826 scaled_positions=scaled_positions,
827 cell=cell)
828 break
830 L += 1
832 # Eigenmodes and -vectors
833 if frequency_factor is None:
834 Kayser_to_eV = 1E2 * 2 * np.pi * units['hbar'] * units['c']
835 # N.B. "fixed default" unit for frequencies in .phonon files is "cm-1"
836 # (i.e. the latter is unaffected by the internal unit conversion system of
837 # CASTEP!) set conversion factor to convert therefrom to eV by default for
838 # now
839 frequency_factor = Kayser_to_eV
840 qpoints = []
841 weights = []
842 frequencies = []
843 displacements = []
844 for _ in range(Nq):
845 fields = lines[L].split()
846 qpoints.append([float(x) for x in fields[2:5]])
847 weights.append(float(fields[5]))
849 freqs = []
850 for _ in range(Nb):
851 L += 1
852 fields = lines[L].split()
853 freqs.append(frequency_factor * float(fields[1]))
854 frequencies.append(np.array(freqs))
856 # skip the two Phonon Eigenvectors header lines
857 L += 2
859 # generate a list of displacements with a structure that is identical to
860 # what is stored internally in the Vibrations class (see in
861 # ase.vibrations.Vibrations.modes):
862 # np.array(displacements).shape == (Nb,3*N)
864 disps = []
865 for _ in range(Nb):
866 disp_coords = []
867 for _ in range(N):
868 L += 1
869 fields = lines[L].split()
870 disp_x = float(fields[2]) + float(fields[3]) * 1.0j
871 disp_y = float(fields[4]) + float(fields[5]) * 1.0j
872 disp_z = float(fields[6]) + float(fields[7]) * 1.0j
873 disp_coords.extend([disp_x, disp_y, disp_z])
874 disps.append(np.array(disp_coords))
875 displacements.append(np.array(disps))
877 L += 1
879 if read_vib_data:
880 if gamma_only:
881 vibdata = [frequencies[0], displacements[0]]
882 else:
883 vibdata = [qpoints, weights, frequencies, displacements]
884 return vibdata, atoms
885 else:
886 return atoms
889# Routines that only the calculator requires
891def read_param(filename='', calc=None, fd=None, get_interface_options=False):
892 if fd is None:
893 if filename == '':
894 raise ValueError('One between filename and fd must be provided')
895 fd = open(filename)
896 elif filename:
897 warnings.warn('Filestream used to read param, file name will be '
898 'ignored')
900 # If necessary, get the interface options
901 if get_interface_options:
902 int_opts = {}
903 optre = re.compile(r'# ASE_INTERFACE ([^\s]+) : ([^\s]+)')
905 lines = fd.readlines()
906 fd.seek(0)
908 for line in lines:
909 m = optre.search(line)
910 if m:
911 int_opts[m.groups()[0]] = m.groups()[1]
913 data = read_freeform(fd)
915 if calc is None:
916 from ase.calculators.castep import Castep
917 calc = Castep(check_castep_version=False, keyword_tolerance=2)
919 for kw, (val, otype) in data.items():
920 if otype == 'block':
921 val = val.split('\n') # Avoids a bug for one-line blocks
922 calc.param.__setattr__(kw, val)
924 if not get_interface_options:
925 return calc
926 else:
927 return calc, int_opts
930def write_param(filename, param, check_checkfile=False,
931 force_write=False,
932 interface_options=None):
933 """Writes a CastepParam object to a CASTEP .param file
935 Parameters
936 ----------
937 filename: the location of the file to write to. If it
938 exists it will be overwritten without warning. If it
939 doesn't it will be created.
940 param: a CastepParam instance
941 check_checkfile : if set to True, write_param will
942 only write continuation or reuse statement
943 if a restart file exists in the same directory
944 """
945 if os.path.isfile(filename) and not force_write:
946 warnings.warn('ase.io.castep.write_param: Set optional argument '
947 'force_write=True to overwrite %s.' % filename)
948 return False
950 out = paropen(filename, 'w')
951 out.write('#######################################################\n')
952 out.write(f'#CASTEP param file: {filename}\n')
953 out.write('#Created using the Atomic Simulation Environment (ASE)#\n')
954 if interface_options is not None:
955 out.write('# Internal settings of the calculator\n')
956 out.write('# This can be switched off by settings\n')
957 out.write('# calc._export_settings = False\n')
958 out.write('# If stated, this will be automatically processed\n')
959 out.write('# by ase.io.castep.read_seed()\n')
960 for option, value in sorted(interface_options.items()):
961 out.write(f'# ASE_INTERFACE {option} : {value}\n')
962 out.write('#######################################################\n\n')
964 if check_checkfile:
965 param = deepcopy(param) # To avoid modifying the parent one
966 for checktype in ['continuation', 'reuse']:
967 opt = getattr(param, checktype)
968 if opt and opt.value:
969 fname = opt.value
970 if fname == 'default':
971 fname = os.path.splitext(filename)[0] + '.check'
972 if not (os.path.exists(fname) or
973 # CASTEP also understands relative path names, hence
974 # also check relative to the param file directory
975 os.path.exists(
976 os.path.join(os.path.dirname(filename),
977 opt.value))):
978 opt.clear()
980 write_freeform(out, param)
982 out.close()
985def read_seed(seed, new_seed=None, ignore_internal_keys=False):
986 """A wrapper around the CASTEP Calculator in conjunction with
987 read_cell and read_param. Basically this can be used to reuse
988 a previous calculation which results in a triple of
989 cell/param/castep file. The label of the calculation if pre-
990 fixed with `copy_of_` and everything else will be recycled as
991 much as possible from the addressed calculation.
993 Please note that this routine will return an atoms ordering as specified
994 in the cell file! It will thus undo the potential reordering internally
995 done by castep.
996 """
998 directory = os.path.abspath(os.path.dirname(seed))
999 seed = os.path.basename(seed)
1001 paramfile = os.path.join(directory, f'{seed}.param')
1002 cellfile = os.path.join(directory, f'{seed}.cell')
1003 castepfile = os.path.join(directory, f'{seed}.castep')
1004 checkfile = os.path.join(directory, f'{seed}.check')
1006 atoms = read_castep_cell(cellfile)
1007 atoms.calc._directory = directory
1008 atoms.calc._rename_existing_dir = False
1009 atoms.calc._castep_pp_path = directory
1010 atoms.calc.merge_param(paramfile,
1011 ignore_internal_keys=ignore_internal_keys)
1012 if new_seed is None:
1013 atoms.calc._label = f'copy_of_{seed}'
1014 else:
1015 atoms.calc._label = str(new_seed)
1016 if os.path.isfile(castepfile):
1017 # _set_atoms needs to be True here
1018 # but we set it right back to False
1019 # atoms.calc._set_atoms = False
1020 # BUGFIX: I do not see a reason to do that!
1021 atoms.calc.read(castepfile)
1022 # atoms.calc._set_atoms = False
1024 # if here is a check file, we also want to re-use this information
1025 if os.path.isfile(checkfile):
1026 atoms.calc._check_file = os.path.basename(checkfile)
1028 # sync the top-level object with the
1029 # one attached to the calculator
1030 atoms = atoms.calc.atoms
1031 else:
1032 # There are cases where we only want to restore a calculator/atoms
1033 # setting without a castep file...
1034 # No print statement required in these cases
1035 warnings.warn(
1036 'Corresponding *.castep file not found. '
1037 'Atoms object will be restored from *.cell and *.param only.')
1038 atoms.calc.push_oldstate()
1040 return atoms
1043@reader
1044def read_bands(fd, units=units_CODATA2002):
1045 """Read Castep.bands file to kpoints, weights and eigenvalues.
1047 Parameters
1048 ----------
1049 fd : str | io.TextIOBase
1050 Path to the `.bands` file or file descriptor for open `.bands` file.
1051 units : dict
1052 Conversion factors for atomic units.
1054 Returns
1055 -------
1056 kpts : np.ndarray
1057 1d NumPy array for k-point coordinates.
1058 weights : np.ndarray
1059 1d NumPy array for k-point weights.
1060 eigenvalues : np.ndarray
1061 NumPy array for eigenvalues with shape (spin, kpts, nbands).
1062 efermi : float
1063 Fermi energy.
1065 """
1066 Hartree = units['Eh']
1068 nkpts = int(fd.readline().split()[-1])
1069 nspin = int(fd.readline().split()[-1])
1070 _ = float(fd.readline().split()[-1])
1071 nbands = int(fd.readline().split()[-1])
1072 efermi = float(fd.readline().split()[-1])
1074 kpts = np.zeros((nkpts, 3))
1075 weights = np.zeros(nkpts)
1076 eigenvalues = np.zeros((nspin, nkpts, nbands))
1078 # Skip unit cell
1079 for _ in range(4):
1080 fd.readline()
1082 def _kptline_to_i_k_wt(line: str) -> tuple[int, list[float], float]:
1083 split_line = line.split()
1084 i_kpt = int(split_line[1]) - 1
1085 kpt = list(map(float, split_line[2:5]))
1086 wt = float(split_line[5])
1087 return i_kpt, kpt, wt
1089 # CASTEP often writes these out-of-order, so check index and write directly
1090 # to the correct row
1091 for _ in range(nkpts):
1092 i_kpt, kpt, wt = _kptline_to_i_k_wt(fd.readline())
1093 kpts[i_kpt, :], weights[i_kpt] = kpt, wt
1094 for spin in range(nspin):
1095 fd.readline() # Skip 'Spin component N' line
1096 eigenvalues[spin, i_kpt, :] = [float(fd.readline())
1097 for _ in range(nbands)]
1099 return (kpts, weights, eigenvalues * Hartree, efermi * Hartree)
1102def _set_spacegroup_info(find_spg: bool, atoms: ase.Atoms) -> None:
1103 """If requested, get spacegroup with spglib and attach to atoms.info"""
1105 if not find_spg:
1106 return
1108 try:
1109 import spglib
1110 except ImportError:
1111 warnings.warn('spglib not found installed on this system - '
1112 'automatic spacegroup detection is not possible')
1113 return
1115 symmd = spglib_new_errorhandling(
1116 spglib.get_symmetry_dataset)(atoms_to_spglib_cell(atoms))
1117 atoms_spg = Spacegroup(int(symmd['number']))
1118 atoms.info['spacegroup'] = atoms_spg