Coverage for /builds/debichem-team/python-ase/ase/io/castep/__init__.py: 50.22%
691 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""This module defines I/O routines with CASTEP files.
2The key idea is that all function accept or return atoms objects.
3CASTEP specific parameters will be returned through the <atoms>.calc
4attribute.
5"""
6import os
7import re
8import warnings
9from copy import deepcopy
10from typing import List, Tuple
12import numpy as np
14import ase
16# independent unit management included here:
17# When high accuracy is required, this allows to easily pin down
18# unit conversion factors from different "unit definition systems"
19# (CODATA1986 for ase-3.6.0.2515 vs CODATA2002 for CASTEP 5.01).
20#
21# ase.units in in ase-3.6.0.2515 is based on CODATA1986
22import ase.units
23from ase.constraints import FixAtoms, FixCartesian, FixedLine, FixedPlane
24from ase.geometry.cell import cellpar_to_cell
25from ase.io.castep.castep_reader import read_castep_castep
26from ase.parallel import paropen
27from ase.spacegroup import Spacegroup
28from ase.utils import atoms_to_spglib_cell, reader, writer
30units_ase = {
31 'hbar': ase.units._hbar * ase.units.J,
32 'Eh': ase.units.Hartree,
33 'kB': ase.units.kB,
34 'a0': ase.units.Bohr,
35 't0': ase.units._hbar * ase.units.J / ase.units.Hartree,
36 'c': ase.units._c,
37 'me': ase.units._me / ase.units._amu,
38 'Pascal': 1.0 / ase.units.Pascal}
40# CODATA1986 (included herein for the sake of completeness)
41# taken from
42# http://physics.nist.gov/cuu/Archive/1986RMP.pdf
43units_CODATA1986 = {
44 'hbar': 6.5821220E-16, # eVs
45 'Eh': 27.2113961, # eV
46 'kB': 8.617385E-5, # eV/K
47 'a0': 0.529177249, # A
48 'c': 299792458, # m/s
49 'e': 1.60217733E-19, # C
50 'me': 5.485799110E-4} # u
52# CODATA2002: default in CASTEP 5.01
53# (-> check in more recent CASTEP in case of numerical discrepancies?!)
54# taken from
55# http://physics.nist.gov/cuu/Document/all_2002.pdf
56units_CODATA2002 = {
57 'hbar': 6.58211915E-16, # eVs
58 'Eh': 27.2113845, # eV
59 'kB': 8.617343E-5, # eV/K
60 'a0': 0.5291772108, # A
61 'c': 299792458, # m/s
62 'e': 1.60217653E-19, # C
63 'me': 5.4857990945E-4} # u
65# (common) derived entries
66for d in (units_CODATA1986, units_CODATA2002):
67 d['t0'] = d['hbar'] / d['Eh'] # s
68 d['Pascal'] = d['e'] * 1E30 # Pa
71__all__ = [
72 # routines for the generic io function
73 'read_castep_castep',
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')
121@writer
122def write_castep_cell(fd, atoms, positions_frac=False, force_write=False,
123 precision=6, magnetic_moments=None,
124 castep_cell=None):
125 """
126 This CASTEP export function write minimal information to
127 a .cell file. If the atoms object is a trajectory, it will
128 take the last image.
130 Note that function has been altered in order to require a filedescriptor
131 rather than a filename. This allows to use the more generic write()
132 function from formats.py
134 Note that the "force_write" keywords has no effect currently.
136 Arguments:
138 positions_frac: boolean. If true, positions are printed as fractional
139 rather than absolute. Default is false.
140 castep_cell: if provided, overrides the existing CastepCell object in
141 the Atoms calculator
142 precision: number of digits to which lattice and positions are printed
143 magnetic_moments: if None, no SPIN values are initialised.
144 If 'initial', the values from
145 get_initial_magnetic_moments() are used.
146 If 'calculated', the values from
147 get_magnetic_moments() are used.
148 If an array of the same length as the atoms object,
149 its contents will be used as magnetic moments.
150 """
152 if isinstance(atoms, list):
153 if len(atoms) > 1:
154 atoms = atoms[-1]
156 # Header
157 fd.write('# written by ASE\n\n')
159 # To write this we simply use the existing Castep calculator, or create
160 # one
161 from ase.calculators.castep import Castep, CastepCell
163 try:
164 has_cell = isinstance(atoms.calc.cell, CastepCell)
165 except AttributeError:
166 has_cell = False
168 if has_cell:
169 cell = deepcopy(atoms.calc.cell)
170 else:
171 cell = Castep(keyword_tolerance=2).cell
173 # Write lattice
174 fformat = f'%{precision + 3}.{precision}f'
175 cell_block_format = ' '.join([fformat] * 3)
176 cell.lattice_cart = [cell_block_format % tuple(line)
177 for line in atoms.get_cell()]
179 if positions_frac:
180 pos_keyword = 'positions_frac'
181 positions = atoms.get_scaled_positions()
182 else:
183 pos_keyword = 'positions_abs'
184 positions = atoms.get_positions()
186 if atoms.has('castep_custom_species'):
187 elems = atoms.get_array('castep_custom_species')
188 else:
189 elems = atoms.get_chemical_symbols()
190 if atoms.has('masses'):
192 from ase.data import atomic_masses
193 masses = atoms.get_array('masses')
194 custom_masses = {}
196 for i, species in enumerate(elems):
197 custom_mass = masses[i]
199 # build record of different masses for each species
200 if species not in custom_masses:
202 # build dictionary of positions of all species with
203 # same name and mass value ideally there should only
204 # be one mass per species
205 custom_masses[species] = {custom_mass: [i]}
207 # if multiple masses found for a species
208 elif custom_mass not in custom_masses[species].keys():
210 # if custom species were already manually defined raise an error
211 if atoms.has('castep_custom_species'):
212 raise ValueError(
213 "Could not write custom mass block for {0}. \n"
214 "Custom mass was set ({1}), but an inconsistent set of "
215 "castep_custom_species already defines "
216 "({2}) for {0}. \n"
217 "If using both features, ensure that "
218 "each species type in "
219 "atoms.arrays['castep_custom_species'] "
220 "has consistent mass values and that each atom "
221 "with non-standard "
222 "mass belongs to a custom species type."
223 "".format(
224 species, custom_mass, list(
225 custom_masses[species].keys())[0]))
227 # append mass to create custom species later
228 else:
229 custom_masses[species][custom_mass] = [i]
230 else:
231 custom_masses[species][custom_mass].append(i)
233 # create species_mass block
234 mass_block = []
236 for el, mass_dict in custom_masses.items():
238 # ignore mass record that match defaults
239 default = mass_dict.pop(atomic_masses[atoms.get_array(
240 'numbers')[list(elems).index(el)]], None)
241 if mass_dict:
242 # no custom species need to be created
243 if len(mass_dict) == 1 and not default:
244 mass_block.append('{} {}'.format(
245 el, list(mass_dict.keys())[0]))
246 # for each custom mass, create new species and change names to
247 # match in 'elems' list
248 else:
249 warnings.warn(
250 'Custom mass specified for '
251 'standard species {}, creating custom species'
252 .format(el))
254 for i, vals in enumerate(mass_dict.items()):
255 mass_val, idxs = vals
256 custom_species_name = f"{el}:{i}"
257 warnings.warn(
258 'Creating custom species {} with mass {}'.format(
259 custom_species_name, str(mass_dict)))
260 for idx in idxs:
261 elems[idx] = custom_species_name
262 mass_block.append('{} {}'.format(
263 custom_species_name, mass_val))
265 cell.species_mass = mass_block
267 if atoms.has('castep_labels'):
268 labels = atoms.get_array('castep_labels')
269 else:
270 labels = ['NULL'] * len(elems)
272 if str(magnetic_moments).lower() == 'initial':
273 magmoms = atoms.get_initial_magnetic_moments()
274 elif str(magnetic_moments).lower() == 'calculated':
275 magmoms = atoms.get_magnetic_moments()
276 elif np.array(magnetic_moments).shape == (len(elems),):
277 magmoms = np.array(magnetic_moments)
278 else:
279 magmoms = [0] * len(elems)
281 pos_block = []
282 pos_block_format = '%s ' + cell_block_format
284 for i, el in enumerate(elems):
285 xyz = positions[i]
286 line = pos_block_format % tuple([el] + list(xyz))
287 # ADD other keywords if necessary
288 if magmoms[i] != 0:
289 line += f' SPIN={magmoms[i]} '
290 if labels[i].strip() not in ('NULL', ''):
291 line += f' LABEL={labels[i]} '
292 pos_block.append(line)
294 setattr(cell, pos_keyword, pos_block)
296 constr_block = _make_block_ionic_constraints(atoms)
297 if constr_block:
298 cell.ionic_constraints = constr_block
300 write_freeform(fd, cell)
303def _make_block_ionic_constraints(atoms: ase.Atoms) -> List[str]:
304 constr_block: List[str] = []
305 species_indices = atoms.symbols.species_indices()
306 for constr in atoms.constraints:
307 if not _is_constraint_valid(constr, len(atoms)):
308 continue
309 for i in constr.index:
310 symbol = atoms.get_chemical_symbols()[i]
311 nis = species_indices[i] + 1
312 if isinstance(constr, FixAtoms):
313 for j in range(3): # constraint for all three directions
314 ic = len(constr_block) + 1
315 line = f'{ic:6d} {symbol:3s} {nis:3d} '
316 line += ['1 0 0', '0 1 0', '0 0 1'][j]
317 constr_block.append(line)
318 elif isinstance(constr, FixCartesian):
319 for j, m in enumerate(constr.mask):
320 if m == 0: # not constrained
321 continue
322 ic = len(constr_block) + 1
323 line = f'{ic:6d} {symbol:3s} {nis:3d} '
324 line += ['1 0 0', '0 1 0', '0 0 1'][j]
325 constr_block.append(line)
326 elif isinstance(constr, FixedPlane):
327 ic = len(constr_block) + 1
328 line = f'{ic:6d} {symbol:3s} {nis:3d} '
329 line += ' '.join([str(d) for d in constr.dir])
330 constr_block.append(line)
331 elif isinstance(constr, FixedLine):
332 for direction in _calc_normal_vectors(constr):
333 ic = len(constr_block) + 1
334 line = f'{ic:6d} {symbol:3s} {nis:3d} '
335 line += ' '.join(str(_) for _ in direction)
336 constr_block.append(line)
337 return constr_block
340def _is_constraint_valid(constraint, natoms: int) -> bool:
341 supported_constraints = (FixAtoms, FixedPlane, FixedLine, FixCartesian)
342 if not isinstance(constraint, supported_constraints):
343 warnings.warn(f'{constraint} is not supported by ASE CASTEP, skipped')
344 return False
345 if any(i < 0 or i >= natoms for i in constraint.index):
346 warnings.warn(f'{constraint} contains invalid indices, skipped')
347 return False
348 return True
351def _calc_normal_vectors(constr: FixedLine) -> Tuple[np.ndarray, np.ndarray]:
352 direction = constr.dir
354 i2, i1 = np.argsort(np.abs(direction))[1:]
355 v1 = direction[i1]
356 v2 = direction[i2]
357 n1 = np.zeros(3)
358 n1[i2] = v1
359 n1[i1] = -v2
360 n1 = n1 / np.linalg.norm(n1)
362 n2 = np.cross(direction, n1)
363 n2 = n2 / np.linalg.norm(n2)
365 return n1, n2
368def read_freeform(fd):
369 """
370 Read a CASTEP freeform file (the basic format of .cell and .param files)
371 and return keyword-value pairs as a dict (values are strings for single
372 keywords and lists of strings for blocks).
373 """
375 from ase.io.castep.castep_input_file import CastepInputFile
377 inputobj = CastepInputFile(keyword_tolerance=2)
379 filelines = fd.readlines()
381 keyw = None
382 read_block = False
383 block_lines = None
385 for i, l in enumerate(filelines):
387 # Strip all comments, aka anything after a hash
388 L = re.split(r'[#!;]', l, 1)[0].strip()
390 if L == '':
391 # Empty line... skip
392 continue
394 lsplit = re.split(r'\s*[:=]*\s+', L, 1)
396 if read_block:
397 if lsplit[0].lower() == '%endblock':
398 if len(lsplit) == 1 or lsplit[1].lower() != keyw:
399 raise ValueError('Out of place end of block at '
400 'line %i in freeform file' % i + 1)
401 else:
402 read_block = False
403 inputobj.__setattr__(keyw, block_lines)
404 else:
405 block_lines += [L]
406 else:
407 # Check the first word
409 # Is it a block?
410 read_block = (lsplit[0].lower() == '%block')
411 if read_block:
412 if len(lsplit) == 1:
413 raise ValueError(('Unrecognizable block at line %i '
414 'in io freeform file') % i + 1)
415 else:
416 keyw = lsplit[1].lower()
417 else:
418 keyw = lsplit[0].lower()
420 # Now save the value
421 if read_block:
422 block_lines = []
423 else:
424 inputobj.__setattr__(keyw, ' '.join(lsplit[1:]))
426 return inputobj.get_attr_dict(types=True)
429@reader
430def read_castep_cell(fd, index=None, calculator_args={}, find_spg=False,
431 units=units_CODATA2002):
432 """Read a .cell file and return an atoms object.
433 Any value found that does not fit the atoms API
434 will be stored in the atoms.calc attribute.
436 By default, the Castep calculator will be tolerant and in the absence of a
437 castep_keywords.json file it will just accept all keywords that aren't
438 automatically parsed.
439 """
441 from ase.calculators.castep import Castep
443 cell_units = { # Units specifiers for CASTEP
444 'bohr': units_CODATA2002['a0'],
445 'ang': 1.0,
446 'm': 1e10,
447 'cm': 1e8,
448 'nm': 10,
449 'pm': 1e-2
450 }
452 calc = Castep(**calculator_args)
454 if calc.cell.castep_version == 0 and calc._kw_tol < 3:
455 # No valid castep_keywords.json was found
456 warnings.warn(
457 'read_cell: Warning - Was not able to validate CASTEP input. '
458 'This may be due to a non-existing '
459 '"castep_keywords.json" '
460 'file or a non-existing CASTEP installation. '
461 'Parsing will go on but keywords will not be '
462 'validated and may cause problems if incorrect during a CASTEP '
463 'run.')
465 celldict = read_freeform(fd)
467 def parse_blockunit(line_tokens, blockname):
468 u = 1.0
469 if len(line_tokens[0]) == 1:
470 usymb = line_tokens[0][0].lower()
471 u = cell_units.get(usymb, 1)
472 if usymb not in cell_units:
473 warnings.warn('read_cell: Warning - ignoring invalid '
474 'unit specifier in %BLOCK {} '
475 '(assuming Angstrom instead)'.format(blockname))
476 line_tokens = line_tokens[1:]
477 return u, line_tokens
479 # Arguments to pass to the Atoms object at the end
480 aargs = {
481 'pbc': True
482 }
484 # Start by looking for the lattice
485 lat_keywords = [w in celldict for w in ('lattice_cart', 'lattice_abc')]
486 if all(lat_keywords):
487 warnings.warn('read_cell: Warning - two lattice blocks present in the'
488 ' same file. LATTICE_ABC will be ignored')
489 elif not any(lat_keywords):
490 raise ValueError('Cell file must contain at least one between '
491 'LATTICE_ABC and LATTICE_CART')
493 if 'lattice_abc' in celldict:
495 lines = celldict.pop('lattice_abc')[0].split('\n')
496 line_tokens = [line.split() for line in lines]
498 u, line_tokens = parse_blockunit(line_tokens, 'lattice_abc')
500 if len(line_tokens) != 2:
501 warnings.warn('read_cell: Warning - ignoring additional '
502 'lines in invalid %BLOCK LATTICE_ABC')
504 abc = [float(p) * u for p in line_tokens[0][:3]]
505 angles = [float(phi) for phi in line_tokens[1][:3]]
507 aargs['cell'] = cellpar_to_cell(abc + angles)
509 if 'lattice_cart' in celldict:
511 lines = celldict.pop('lattice_cart')[0].split('\n')
512 line_tokens = [line.split() for line in lines]
514 u, line_tokens = parse_blockunit(line_tokens, 'lattice_cart')
516 if len(line_tokens) != 3:
517 warnings.warn('read_cell: Warning - ignoring more than '
518 'three lattice vectors in invalid %BLOCK '
519 'LATTICE_CART')
521 aargs['cell'] = [[float(x) * u for x in lt[:3]] for lt in line_tokens]
523 # Now move on to the positions
524 pos_keywords = [w in celldict
525 for w in ('positions_abs', 'positions_frac')]
527 if all(pos_keywords):
528 warnings.warn('read_cell: Warning - two lattice blocks present in the'
529 ' same file. POSITIONS_FRAC will be ignored')
530 del celldict['positions_frac']
531 elif not any(pos_keywords):
532 raise ValueError('Cell file must contain at least one between '
533 'POSITIONS_FRAC and POSITIONS_ABS')
535 aargs['symbols'] = []
536 pos_type = 'positions'
537 pos_block = celldict.pop('positions_abs', [None])[0]
538 if pos_block is None:
539 pos_type = 'scaled_positions'
540 pos_block = celldict.pop('positions_frac', [None])[0]
541 aargs[pos_type] = []
543 lines = pos_block.split('\n')
544 line_tokens = [line.split() for line in lines]
546 if 'scaled' not in pos_type:
547 u, line_tokens = parse_blockunit(line_tokens, 'positions_abs')
548 else:
549 u = 1.0
551 # Here we extract all the possible additional info
552 # These are marked by their type
554 add_info = {
555 'SPIN': (float, 0.0), # (type, default)
556 'MAGMOM': (float, 0.0),
557 'LABEL': (str, 'NULL')
558 }
559 add_info_arrays = {k: [] for k in add_info}
561 def parse_info(raw_info):
563 re_keys = (r'({0})\s*[=:\s]{{1}}\s'
564 r'*([^\s]*)').format('|'.join(add_info.keys()))
565 # Capture all info groups
566 info = re.findall(re_keys, raw_info)
567 info = {g[0]: add_info[g[0]][0](g[1]) for g in info}
568 return info
570 # Array for custom species (a CASTEP special thing)
571 # Usually left unused
572 custom_species = None
574 for tokens in line_tokens:
575 # Now, process the whole 'species' thing
576 spec_custom = tokens[0].split(':', 1)
577 elem = spec_custom[0]
578 if len(spec_custom) > 1 and custom_species is None:
579 # Add it to the custom info!
580 custom_species = list(aargs['symbols'])
581 if custom_species is not None:
582 custom_species.append(tokens[0])
583 aargs['symbols'].append(elem)
584 aargs[pos_type].append([float(p) * u for p in tokens[1:4]])
585 # Now for the additional information
586 info = ' '.join(tokens[4:])
587 info = parse_info(info)
588 for k in add_info:
589 add_info_arrays[k] += [info.get(k, add_info[k][1])]
591 # read in custom species mass
592 if 'species_mass' in celldict:
593 spec_list = custom_species if custom_species else aargs['symbols']
594 aargs['masses'] = [None for _ in spec_list]
595 lines = celldict.pop('species_mass')[0].split('\n')
596 line_tokens = [line.split() for line in lines]
598 if len(line_tokens[0]) == 1:
599 if line_tokens[0][0].lower() not in ('amu', 'u'):
600 raise ValueError(
601 "unit specifier '{}' in %BLOCK SPECIES_MASS "
602 "not recognised".format(
603 line_tokens[0][0].lower()))
604 line_tokens = line_tokens[1:]
606 for tokens in line_tokens:
607 token_pos_list = [i for i, x in enumerate(
608 spec_list) if x == tokens[0]]
609 if len(token_pos_list) == 0:
610 warnings.warn(
611 'read_cell: Warning - ignoring unused '
612 'species mass {} in %BLOCK SPECIES_MASS'.format(
613 tokens[0]))
614 for idx in token_pos_list:
615 aargs['masses'][idx] = tokens[1]
617 # Now on to the species potentials...
618 if 'species_pot' in celldict:
619 lines = celldict.pop('species_pot')[0].split('\n')
620 line_tokens = [line.split() for line in lines]
622 for tokens in line_tokens:
623 if len(tokens) == 1:
624 # It's a library
625 all_spec = (set(custom_species) if custom_species is not None
626 else set(aargs['symbols']))
627 for s in all_spec:
628 calc.cell.species_pot = (s, tokens[0])
629 else:
630 calc.cell.species_pot = tuple(tokens[:2])
632 # Ionic constraints
633 raw_constraints = {}
635 if 'ionic_constraints' in celldict:
636 lines = celldict.pop('ionic_constraints')[0].split('\n')
637 line_tokens = [line.split() for line in lines]
639 for tokens in line_tokens:
640 if len(tokens) != 6:
641 continue
642 _, species, nic, x, y, z = tokens
643 # convert xyz to floats
644 x = float(x)
645 y = float(y)
646 z = float(z)
648 nic = int(nic)
649 if (species, nic) not in raw_constraints:
650 raw_constraints[(species, nic)] = []
651 raw_constraints[(species, nic)].append(np.array(
652 [x, y, z]))
654 # Symmetry operations
655 if 'symmetry_ops' in celldict:
656 lines = celldict.pop('symmetry_ops')[0].split('\n')
657 line_tokens = [line.split() for line in lines]
659 # Read them in blocks of four
660 blocks = np.array(line_tokens).astype(float)
661 if (len(blocks.shape) != 2 or blocks.shape[1] != 3
662 or blocks.shape[0] % 4 != 0):
663 warnings.warn('Warning: could not parse SYMMETRY_OPS'
664 ' block properly, skipping')
665 else:
666 blocks = blocks.reshape((-1, 4, 3))
667 rotations = blocks[:, :3]
668 translations = blocks[:, 3]
670 # Regardless of whether we recognize them, store these
671 calc.cell.symmetry_ops = (rotations, translations)
673 # Anything else that remains, just add it to the cell object:
674 for k, (val, otype) in celldict.items():
675 try:
676 if otype == 'block':
677 val = val.split('\n') # Avoids a bug for one-line blocks
678 calc.cell.__setattr__(k, val)
679 except Exception as e:
680 raise RuntimeError(
681 f'Problem setting calc.cell.{k} = {val}: {e}')
683 # Get the relevant additional info
684 aargs['magmoms'] = np.array(add_info_arrays['SPIN'])
685 # SPIN or MAGMOM are alternative keywords
686 aargs['magmoms'] = np.where(aargs['magmoms'] != 0,
687 aargs['magmoms'],
688 add_info_arrays['MAGMOM'])
689 labels = np.array(add_info_arrays['LABEL'])
691 aargs['calculator'] = calc
693 atoms = ase.Atoms(**aargs)
695 # Spacegroup...
696 if find_spg:
697 # Try importing spglib
698 try:
699 import spglib
700 except ImportError:
701 warnings.warn('spglib not found installed on this system - '
702 'automatic spacegroup detection is not possible')
703 spglib = None
705 if spglib is not None:
706 symmd = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms))
707 atoms_spg = Spacegroup(int(symmd['number']))
708 atoms.info['spacegroup'] = atoms_spg
710 atoms.new_array('castep_labels', labels)
711 if custom_species is not None:
712 atoms.new_array('castep_custom_species', np.array(custom_species))
714 fixed_atoms = []
715 constraints = []
716 index_dict = atoms.symbols.indices()
717 for (species, nic), value in raw_constraints.items():
719 absolute_nr = index_dict[species][nic - 1]
720 if len(value) == 3:
721 # Check if they are linearly independent
722 if np.linalg.det(value) == 0:
723 warnings.warn(
724 'Error: Found linearly dependent constraints attached '
725 'to atoms %s' %
726 (absolute_nr))
727 continue
728 fixed_atoms.append(absolute_nr)
729 elif len(value) == 2:
730 direction = np.cross(value[0], value[1])
731 # Check if they are linearly independent
732 if np.linalg.norm(direction) == 0:
733 warnings.warn(
734 'Error: Found linearly dependent constraints attached '
735 'to atoms %s' %
736 (absolute_nr))
737 continue
738 constraint = FixedLine(indices=absolute_nr, direction=direction)
739 constraints.append(constraint)
740 elif len(value) == 1:
741 direction = np.array(value[0], dtype=float)
742 constraint = FixedPlane(indices=absolute_nr, direction=direction)
743 constraints.append(constraint)
744 else:
745 warnings.warn(
746 f'Error: Found {len(value)} statements attached to atoms '
747 f'{absolute_nr}'
748 )
750 # we need to sort the fixed atoms list in order not to raise an assertion
751 # error in FixAtoms
752 if fixed_atoms:
753 constraints.append(FixAtoms(indices=sorted(fixed_atoms)))
754 if constraints:
755 atoms.set_constraint(constraints)
757 atoms.calc.atoms = atoms
758 atoms.calc.push_oldstate()
760 return atoms
763def read_geom(filename, index=':', units=units_CODATA2002):
764 """
765 Wrapper function for the more generic read() functionality.
767 Note that this is function is intended to maintain backwards-compatibility
768 only. Keyword arguments will be passed to read_castep_geom().
769 """
770 from ase.io import read
771 return read(filename, index=index, format='castep-geom', units=units)
774def read_castep_geom(fd, index=None, units=units_CODATA2002):
775 """Reads a .geom file produced by the CASTEP GeometryOptimization task and
776 returns an atoms object.
777 The information about total free energy and forces of each atom for every
778 relaxation step will be stored for further analysis especially in a
779 single-point calculator.
780 Note that everything in the .geom file is in atomic units, which has
781 been conversed to commonly used unit angstrom(length) and eV (energy).
783 Note that the index argument has no effect as of now.
785 Contribution by Wei-Bing Zhang. Thanks!
787 Routine now accepts a filedescriptor in order to out-source the gz and
788 bz2 handling to formats.py. Note that there is a fallback routine
789 read_geom() that behaves like previous versions did.
790 """
791 from ase.calculators.singlepoint import SinglePointCalculator
793 # fd is closed by embracing read() routine
794 txt = fd.readlines()
796 traj = []
798 Hartree = units['Eh']
799 Bohr = units['a0']
801 # Yeah, we know that...
802 # print('N.B.: Energy in .geom file is not 0K extrapolated.')
803 for i, line in enumerate(txt):
804 if line.find('<-- E') > 0:
805 start_found = True
806 energy = float(line.split()[0]) * Hartree
807 cell = [x.split()[0:3] for x in txt[i + 1:i + 4]]
808 cell = np.array([[float(col) * Bohr for col in row] for row in
809 cell])
810 if line.find('<-- R') > 0 and start_found:
811 start_found = False
812 geom_start = i
813 for i, line in enumerate(txt[geom_start:]):
814 if line.find('<-- F') > 0:
815 geom_stop = i + geom_start
816 break
817 species = [line.split()[0] for line in
818 txt[geom_start:geom_stop]]
819 geom = np.array([[float(col) * Bohr for col in
820 line.split()[2:5]] for line in
821 txt[geom_start:geom_stop]])
822 forces = np.array([[float(col) * Hartree / Bohr for col in
823 line.split()[2:5]] for line in
824 txt[geom_stop:geom_stop
825 + (geom_stop - geom_start)]])
826 image = ase.Atoms(species, geom, cell=cell, pbc=True)
827 image.calc = SinglePointCalculator(
828 atoms=image, energy=energy, forces=forces)
829 traj.append(image)
831 if index is None:
832 return traj
833 else:
834 return traj[index]
837def read_phonon(filename, index=None, read_vib_data=False,
838 gamma_only=True, frequency_factor=None,
839 units=units_CODATA2002):
840 """
841 Wrapper function for the more generic read() functionality.
843 Note that this is function is intended to maintain backwards-compatibility
844 only. For documentation see read_castep_phonon().
845 """
846 from ase.io import read
848 if read_vib_data:
849 full_output = True
850 else:
851 full_output = False
853 return read(filename, index=index, format='castep-phonon',
854 full_output=full_output, read_vib_data=read_vib_data,
855 gamma_only=gamma_only, frequency_factor=frequency_factor,
856 units=units)
859def read_castep_phonon(fd, index=None, read_vib_data=False,
860 gamma_only=True, frequency_factor=None,
861 units=units_CODATA2002):
862 """
863 Reads a .phonon file written by a CASTEP Phonon task and returns an atoms
864 object, as well as the calculated vibrational data if requested.
866 Note that the index argument has no effect as of now.
867 """
869 # fd is closed by embracing read() routine
870 lines = fd.readlines()
872 atoms = None
873 cell = []
874 N = Nb = Nq = 0
875 scaled_positions = []
876 symbols = []
877 masses = []
879 # header
880 L = 0
881 while L < len(lines):
883 line = lines[L]
885 if 'Number of ions' in line:
886 N = int(line.split()[3])
887 elif 'Number of branches' in line:
888 Nb = int(line.split()[3])
889 elif 'Number of wavevectors' in line:
890 Nq = int(line.split()[3])
891 elif 'Unit cell vectors (A)' in line:
892 for _ in range(3):
893 L += 1
894 fields = lines[L].split()
895 cell.append([float(x) for x in fields[0:3]])
896 elif 'Fractional Co-ordinates' in line:
897 for _ in range(N):
898 L += 1
899 fields = lines[L].split()
900 scaled_positions.append([float(x) for x in fields[1:4]])
901 symbols.append(fields[4])
902 masses.append(float(fields[5]))
903 elif 'END header' in line:
904 L += 1
905 atoms = ase.Atoms(symbols=symbols,
906 scaled_positions=scaled_positions,
907 cell=cell)
908 break
910 L += 1
912 # Eigenmodes and -vectors
913 if frequency_factor is None:
914 Kayser_to_eV = 1E2 * 2 * np.pi * units['hbar'] * units['c']
915 # N.B. "fixed default" unit for frequencies in .phonon files is "cm-1"
916 # (i.e. the latter is unaffected by the internal unit conversion system of
917 # CASTEP!) set conversion factor to convert therefrom to eV by default for
918 # now
919 frequency_factor = Kayser_to_eV
920 qpoints = []
921 weights = []
922 frequencies = []
923 displacements = []
924 for _ in range(Nq):
925 fields = lines[L].split()
926 qpoints.append([float(x) for x in fields[2:5]])
927 weights.append(float(fields[5]))
928 freqs = []
929 for _ in range(Nb):
930 L += 1
931 fields = lines[L].split()
932 freqs.append(frequency_factor * float(fields[1]))
933 frequencies.append(np.array(freqs))
935 # skip the two Phonon Eigenvectors header lines
936 L += 2
938 # generate a list of displacements with a structure that is identical to
939 # what is stored internally in the Vibrations class (see in
940 # ase.vibrations.Vibrations.modes):
941 # np.array(displacements).shape == (Nb,3*N)
943 disps = []
944 for _ in range(Nb):
945 disp_coords = []
946 for _ in range(N):
947 L += 1
948 fields = lines[L].split()
949 disp_x = float(fields[2]) + float(fields[3]) * 1.0j
950 disp_y = float(fields[4]) + float(fields[5]) * 1.0j
951 disp_z = float(fields[6]) + float(fields[7]) * 1.0j
952 disp_coords.extend([disp_x, disp_y, disp_z])
953 disps.append(np.array(disp_coords))
954 displacements.append(np.array(disps))
956 if read_vib_data:
957 if gamma_only:
958 vibdata = [frequencies[0], displacements[0]]
959 else:
960 vibdata = [qpoints, weights, frequencies, displacements]
961 return vibdata, atoms
962 else:
963 return atoms
966def read_md(filename, index=None, return_scalars=False,
967 units=units_CODATA2002):
968 """Wrapper function for the more generic read() functionality.
970 Note that this function is intended to maintain backwards-compatibility
971 only. For documentation see read_castep_md()
972 """
973 if return_scalars:
974 full_output = True
975 else:
976 full_output = False
978 from ase.io import read
979 return read(filename, index=index, format='castep-md',
980 full_output=full_output, return_scalars=return_scalars,
981 units=units)
984def read_castep_md(fd, index=None, return_scalars=False,
985 units=units_CODATA2002):
986 """Reads a .md file written by a CASTEP MolecularDynamics task
987 and returns the trajectory stored therein as a list of atoms object.
989 Note that the index argument has no effect as of now."""
991 from ase.calculators.singlepoint import SinglePointCalculator
993 factors = {
994 't': units['t0'] * 1E15, # fs
995 'E': units['Eh'], # eV
996 'T': units['Eh'] / units['kB'],
997 'P': units['Eh'] / units['a0']**3 * units['Pascal'],
998 'h': units['a0'],
999 'hv': units['a0'] / units['t0'],
1000 'S': units['Eh'] / units['a0']**3,
1001 'R': units['a0'],
1002 'V': np.sqrt(units['Eh'] / units['me']),
1003 'F': units['Eh'] / units['a0']}
1005 # fd is closed by embracing read() routine
1006 lines = fd.readlines()
1008 L = 0
1009 while 'END header' not in lines[L]:
1010 L += 1
1011 l_end_header = L
1012 lines = lines[l_end_header + 1:]
1013 times = []
1014 energies = []
1015 temperatures = []
1016 pressures = []
1017 traj = []
1019 # Initialization
1020 time = None
1021 Epot = None
1022 Ekin = None
1023 EH = None
1024 temperature = None
1025 pressure = None
1026 symbols = None
1027 positions = None
1028 cell = None
1029 velocities = None
1030 symbols = []
1031 positions = []
1032 velocities = []
1033 forces = []
1034 cell = np.eye(3)
1035 cell_velocities = []
1036 stress = []
1038 for (L, line) in enumerate(lines):
1039 fields = line.split()
1040 if len(fields) == 0:
1041 if L != 0:
1042 times.append(time)
1043 energies.append([Epot, EH, Ekin])
1044 temperatures.append(temperature)
1045 pressures.append(pressure)
1046 atoms = ase.Atoms(symbols=symbols,
1047 positions=positions,
1048 cell=cell)
1049 atoms.set_velocities(velocities)
1050 if len(stress) == 0:
1051 atoms.calc = SinglePointCalculator(
1052 atoms=atoms, energy=Epot, forces=forces)
1053 else:
1054 atoms.calc = SinglePointCalculator(
1055 atoms=atoms, energy=Epot,
1056 forces=forces, stress=stress)
1057 traj.append(atoms)
1058 symbols = []
1059 positions = []
1060 velocities = []
1061 forces = []
1062 cell = []
1063 cell_velocities = []
1064 stress = []
1065 continue
1066 if len(fields) == 1:
1067 time = factors['t'] * float(fields[0])
1068 continue
1070 if fields[-1] == 'E':
1071 E = [float(x) for x in fields[0:3]]
1072 Epot, EH, Ekin = (factors['E'] * Ei for Ei in E)
1073 continue
1075 if fields[-1] == 'T':
1076 temperature = factors['T'] * float(fields[0])
1077 continue
1079 # only printed in case of variable cell calculation or calculate_stress
1080 # explicitly requested
1081 if fields[-1] == 'P':
1082 pressure = factors['P'] * float(fields[0])
1083 continue
1084 if fields[-1] == 'h':
1085 h = [float(x) for x in fields[0:3]]
1086 cell.append([factors['h'] * hi for hi in h])
1087 continue
1089 # only printed in case of variable cell calculation
1090 if fields[-1] == 'hv':
1091 hv = [float(x) for x in fields[0:3]]
1092 cell_velocities.append([factors['hv'] * hvi for hvi in hv])
1093 continue
1095 # only printed in case of variable cell calculation
1096 if fields[-1] == 'S':
1097 S = [float(x) for x in fields[0:3]]
1098 stress.append([factors['S'] * Si for Si in S])
1099 continue
1100 if fields[-1] == 'R':
1101 symbols.append(fields[0])
1102 R = [float(x) for x in fields[2:5]]
1103 positions.append([factors['R'] * Ri for Ri in R])
1104 continue
1105 if fields[-1] == 'V':
1106 V = [float(x) for x in fields[2:5]]
1107 velocities.append([factors['V'] * Vi for Vi in V])
1108 continue
1109 if fields[-1] == 'F':
1110 F = [float(x) for x in fields[2:5]]
1111 forces.append([factors['F'] * Fi for Fi in F])
1112 continue
1114 if index is None:
1115 pass
1116 else:
1117 traj = traj[index]
1119 if return_scalars:
1120 data = [times, energies, temperatures, pressures]
1121 return data, traj
1122 else:
1123 return traj
1126# Routines that only the calculator requires
1128def read_param(filename='', calc=None, fd=None, get_interface_options=False):
1129 if fd is None:
1130 if filename == '':
1131 raise ValueError('One between filename and fd must be provided')
1132 fd = open(filename)
1133 elif filename:
1134 warnings.warn('Filestream used to read param, file name will be '
1135 'ignored')
1137 # If necessary, get the interface options
1138 if get_interface_options:
1139 int_opts = {}
1140 optre = re.compile(r'# ASE_INTERFACE ([^\s]+) : ([^\s]+)')
1142 lines = fd.readlines()
1143 fd.seek(0)
1145 for line in lines:
1146 m = optre.search(line)
1147 if m:
1148 int_opts[m.groups()[0]] = m.groups()[1]
1150 data = read_freeform(fd)
1152 if calc is None:
1153 from ase.calculators.castep import Castep
1154 calc = Castep(check_castep_version=False, keyword_tolerance=2)
1156 for kw, (val, otype) in data.items():
1157 if otype == 'block':
1158 val = val.split('\n') # Avoids a bug for one-line blocks
1159 calc.param.__setattr__(kw, val)
1161 if not get_interface_options:
1162 return calc
1163 else:
1164 return calc, int_opts
1167def write_param(filename, param, check_checkfile=False,
1168 force_write=False,
1169 interface_options=None):
1170 """Writes a CastepParam object to a CASTEP .param file
1172 Parameters:
1173 filename: the location of the file to write to. If it
1174 exists it will be overwritten without warning. If it
1175 doesn't it will be created.
1176 param: a CastepParam instance
1177 check_checkfile : if set to True, write_param will
1178 only write continuation or reuse statement
1179 if a restart file exists in the same directory
1180 """
1181 if os.path.isfile(filename) and not force_write:
1182 warnings.warn('ase.io.castep.write_param: Set optional argument '
1183 'force_write=True to overwrite %s.' % filename)
1184 return False
1186 out = paropen(filename, 'w')
1187 out.write('#######################################################\n')
1188 out.write(f'#CASTEP param file: {filename}\n')
1189 out.write('#Created using the Atomic Simulation Environment (ASE)#\n')
1190 if interface_options is not None:
1191 out.write('# Internal settings of the calculator\n')
1192 out.write('# This can be switched off by settings\n')
1193 out.write('# calc._export_settings = False\n')
1194 out.write('# If stated, this will be automatically processed\n')
1195 out.write('# by ase.io.castep.read_seed()\n')
1196 for option, value in sorted(interface_options.items()):
1197 out.write(f'# ASE_INTERFACE {option} : {value}\n')
1198 out.write('#######################################################\n\n')
1200 if check_checkfile:
1201 param = deepcopy(param) # To avoid modifying the parent one
1202 for checktype in ['continuation', 'reuse']:
1203 opt = getattr(param, checktype)
1204 if opt and opt.value:
1205 fname = opt.value
1206 if fname == 'default':
1207 fname = os.path.splitext(filename)[0] + '.check'
1208 if not (os.path.exists(fname) or
1209 # CASTEP also understands relative path names, hence
1210 # also check relative to the param file directory
1211 os.path.exists(
1212 os.path.join(os.path.dirname(filename),
1213 opt.value))):
1214 opt.clear()
1216 write_freeform(out, param)
1218 out.close()
1221def read_seed(seed, new_seed=None, ignore_internal_keys=False):
1222 """A wrapper around the CASTEP Calculator in conjunction with
1223 read_cell and read_param. Basically this can be used to reuse
1224 a previous calculation which results in a triple of
1225 cell/param/castep file. The label of the calculation if pre-
1226 fixed with `copy_of_` and everything else will be recycled as
1227 much as possible from the addressed calculation.
1229 Please note that this routine will return an atoms ordering as specified
1230 in the cell file! It will thus undo the potential reordering internally
1231 done by castep.
1232 """
1234 directory = os.path.abspath(os.path.dirname(seed))
1235 seed = os.path.basename(seed)
1237 paramfile = os.path.join(directory, f'{seed}.param')
1238 cellfile = os.path.join(directory, f'{seed}.cell')
1239 castepfile = os.path.join(directory, f'{seed}.castep')
1240 checkfile = os.path.join(directory, f'{seed}.check')
1242 atoms = read_castep_cell(cellfile)
1243 atoms.calc._directory = directory
1244 atoms.calc._rename_existing_dir = False
1245 atoms.calc._castep_pp_path = directory
1246 atoms.calc.merge_param(paramfile,
1247 ignore_internal_keys=ignore_internal_keys)
1248 if new_seed is None:
1249 atoms.calc._label = f'copy_of_{seed}'
1250 else:
1251 atoms.calc._label = str(new_seed)
1252 if os.path.isfile(castepfile):
1253 # _set_atoms needs to be True here
1254 # but we set it right back to False
1255 # atoms.calc._set_atoms = False
1256 # BUGFIX: I do not see a reason to do that!
1257 atoms.calc.read(castepfile)
1258 # atoms.calc._set_atoms = False
1260 # if here is a check file, we also want to re-use this information
1261 if os.path.isfile(checkfile):
1262 atoms.calc._check_file = os.path.basename(checkfile)
1264 # sync the top-level object with the
1265 # one attached to the calculator
1266 atoms = atoms.calc.atoms
1267 else:
1268 # There are cases where we only want to restore a calculator/atoms
1269 # setting without a castep file...
1270 # No print statement required in these cases
1271 warnings.warn(
1272 'Corresponding *.castep file not found. '
1273 'Atoms object will be restored from *.cell and *.param only.')
1274 atoms.calc.push_oldstate()
1276 return atoms
1279@reader
1280def read_bands(fd, units=units_CODATA2002):
1281 """Read Castep.bands file to kpoints, weights and eigenvalues.
1283 Parameters
1284 ----------
1285 fd : str | io.TextIOBase
1286 Path to the `.bands` file or file descriptor for open `.bands` file.
1287 units : dict
1288 Conversion factors for atomic units.
1290 Returns
1291 -------
1292 kpts : np.ndarray
1293 1d NumPy array for k-point coordinates.
1294 weights : np.ndarray
1295 1d NumPy array for k-point weights.
1296 eigenvalues : np.ndarray
1297 NumPy array for eigenvalues with shape (spin, kpts, nbands).
1298 efermi : float
1299 Fermi energy.
1301 """
1302 Hartree = units['Eh']
1304 nkpts = int(fd.readline().split()[-1])
1305 nspin = int(fd.readline().split()[-1])
1306 _ = float(fd.readline().split()[-1])
1307 nbands = int(fd.readline().split()[-1])
1308 efermi = float(fd.readline().split()[-1])
1310 kpts = np.zeros((nkpts, 3))
1311 weights = np.zeros(nkpts)
1312 eigenvalues = np.zeros((nspin, nkpts, nbands))
1314 # Skip unit cell
1315 for _ in range(4):
1316 fd.readline()
1318 def _kptline_to_i_k_wt(line: str) -> Tuple[int, List[float], float]:
1319 split_line = line.split()
1320 i_kpt = int(split_line[1]) - 1
1321 kpt = list(map(float, split_line[2:5]))
1322 wt = float(split_line[5])
1323 return i_kpt, kpt, wt
1325 # CASTEP often writes these out-of-order, so check index and write directly
1326 # to the correct row
1327 for _ in range(nkpts):
1328 i_kpt, kpt, wt = _kptline_to_i_k_wt(fd.readline())
1329 kpts[i_kpt, :], weights[i_kpt] = kpt, wt
1330 for spin in range(nspin):
1331 fd.readline() # Skip 'Spin component N' line
1332 eigenvalues[spin, i_kpt, :] = [float(fd.readline())
1333 for _ in range(nbands)]
1335 return (kpts, weights, eigenvalues * Hartree, efermi * Hartree)