Coverage for /builds/debichem-team/python-ase/ase/io/castep.py : 29.36%
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()
210 if atoms.has('castep_labels'):
211 labels = atoms.get_array('castep_labels')
212 else:
213 labels = ['NULL'] * len(elems)
215 if str(magnetic_moments).lower() == 'initial':
216 magmoms = atoms.get_initial_magnetic_moments()
217 elif str(magnetic_moments).lower() == 'calculated':
218 magmoms = atoms.get_magnetic_moments()
219 elif np.array(magnetic_moments).shape == (len(elems),):
220 magmoms = np.array(magnetic_moments)
221 else:
222 magmoms = [0] * len(elems)
224 pos_block = []
225 pos_block_format = '%s ' + cell_block_format
227 for i, el in enumerate(elems):
228 xyz = positions[i]
229 line = pos_block_format % tuple([el] + list(xyz))
230 # ADD other keywords if necessary
231 if magmoms[i] != 0:
232 line += ' SPIN={0} '.format(magmoms[i])
233 if labels[i].strip() not in ('NULL', ''):
234 line += ' LABEL={0} '.format(labels[i])
235 pos_block.append(line)
237 setattr(cell, pos_keyword, pos_block)
239 constraints = atoms.constraints
240 if len(constraints):
241 _supported_constraints = (FixAtoms, FixedPlane, FixedLine,
242 FixCartesian)
244 constr_block = []
246 for constr in constraints:
247 if not isinstance(constr, _supported_constraints):
248 warnings.warn('Warning: you have constraints in your atoms, that are '
249 'not supported by the CASTEP ase interface')
250 break
251 if isinstance(constr, FixAtoms):
252 for i in constr.index:
254 try:
255 symbol = atoms.get_chemical_symbols()[i]
256 nis = atoms.calc._get_number_in_species(i)
257 except KeyError:
258 raise UserWarning('Unrecognized index in'
259 + ' constraint %s' % constr)
260 for j in range(3):
261 L = '%6d %3s %3d ' % (len(constr_block) + 1,
262 symbol,
263 nis)
264 L += ['1 0 0', '0 1 0', '0 0 1'][j]
265 constr_block += [L]
267 elif isinstance(constr, FixCartesian):
268 n = constr.a
269 symbol = atoms.get_chemical_symbols()[n]
270 nis = atoms.calc._get_number_in_species(n)
272 for i, m in enumerate(constr.mask):
273 if m == 1:
274 continue
275 L = '%6d %3s %3d ' % (len(constr_block) + 1, symbol, nis)
276 L += ' '.join(['1' if j == i else '0' for j in range(3)])
277 constr_block += [L]
279 elif isinstance(constr, FixedPlane):
280 n = constr.a
281 symbol = atoms.get_chemical_symbols()[n]
282 nis = atoms.calc._get_number_in_species(n)
284 L = '%6d %3s %3d ' % (len(constr_block) + 1, symbol, nis)
285 L += ' '.join([str(d) for d in constr.dir])
286 constr_block += [L]
288 elif isinstance(constr, FixedLine):
289 n = constr.a
290 symbol = atoms.get_chemical_symbols()[n]
291 nis = atoms.calc._get_number_in_species(n)
293 direction = constr.dir
294 ((i1, v1), (i2, v2)) = sorted(enumerate(direction),
295 key=lambda x: abs(x[1]),
296 reverse=True)[:2]
297 n1 = np.zeros(3)
298 n1[i2] = v1
299 n1[i1] = -v2
300 n1 = n1 / np.linalg.norm(n1)
302 n2 = np.cross(direction, n1)
304 l1 = '%6d %3s %3d %f %f %f' % (len(constr_block) + 1,
305 symbol, nis,
306 n1[0], n1[1], n1[2])
307 l2 = '%6d %3s %3d %f %f %f' % (len(constr_block) + 2,
308 symbol, nis,
309 n2[0], n2[1], n2[2])
311 constr_block += [l1, l2]
313 cell.ionic_constraints = constr_block
315 write_freeform(fd, cell)
317 return True
320def read_freeform(fd):
321 """
322 Read a CASTEP freeform file (the basic format of .cell and .param files)
323 and return keyword-value pairs as a dict (values are strings for single
324 keywords and lists of strings for blocks).
325 """
327 from ase.calculators.castep import CastepInputFile
329 inputobj = CastepInputFile(keyword_tolerance=2)
331 filelines = fd.readlines()
333 keyw = None
334 read_block = False
335 block_lines = None
337 for i, l in enumerate(filelines):
339 # Strip all comments, aka anything after a hash
340 L = re.split(r'[#!;]', l, 1)[0].strip()
342 if L == '':
343 # Empty line... skip
344 continue
346 lsplit = re.split(r'\s*[:=]*\s+', L, 1)
348 if read_block:
349 if lsplit[0].lower() == '%endblock':
350 if len(lsplit) == 1 or lsplit[1].lower() != keyw:
351 raise ValueError('Out of place end of block at '
352 'line %i in freeform file' % i + 1)
353 else:
354 read_block = False
355 inputobj.__setattr__(keyw, block_lines)
356 else:
357 block_lines += [L]
358 else:
359 # Check the first word
361 # Is it a block?
362 read_block = (lsplit[0].lower() == '%block')
363 if read_block:
364 if len(lsplit) == 1:
365 raise ValueError(('Unrecognizable block at line %i '
366 'in io freeform file') % i + 1)
367 else:
368 keyw = lsplit[1].lower()
369 else:
370 keyw = lsplit[0].lower()
372 # Now save the value
373 if read_block:
374 block_lines = []
375 else:
376 inputobj.__setattr__(keyw, ' '.join(lsplit[1:]))
378 return inputobj.get_attr_dict(types=True)
381def read_cell(filename, index=None):
382 """
383 Wrapper function for the more generic read() functionality.
385 Note that this is function is intended to maintain backwards-compatibility
386 only.
387 """
388 from ase.io import read
389 return read(filename, index=index, format='castep-cell')
392def read_castep_cell(fd, index=None, calculator_args={}, find_spg=False,
393 units=units_CODATA2002):
394 """Read a .cell file and return an atoms object.
395 Any value found that does not fit the atoms API
396 will be stored in the atoms.calc attribute.
398 By default, the Castep calculator will be tolerant and in the absence of a
399 castep_keywords.json file it will just accept all keywords that aren't
400 automatically parsed.
401 """
403 from ase.calculators.castep import Castep
405 cell_units = { # Units specifiers for CASTEP
406 'bohr': units_CODATA2002['a0'],
407 'ang': 1.0,
408 'm': 1e10,
409 'cm': 1e8,
410 'nm': 10,
411 'pm': 1e-2
412 }
414 calc = Castep(**calculator_args)
416 if calc.cell.castep_version == 0 and calc._kw_tol < 3:
417 # No valid castep_keywords.json was found
418 warnings.warn('read_cell: Warning - Was not able to validate CASTEP input. '
419 'This may be due to a non-existing '
420 '"castep_keywords.json" '
421 'file or a non-existing CASTEP installation. '
422 'Parsing will go on but keywords will not be '
423 'validated and may cause problems if incorrect during a CASTEP '
424 'run.')
426 celldict = read_freeform(fd)
428 def parse_blockunit(line_tokens, blockname):
429 u = 1.0
430 if len(line_tokens[0]) == 1:
431 usymb = line_tokens[0][0].lower()
432 u = cell_units.get(usymb, 1)
433 if usymb not in cell_units:
434 warnings.warn('read_cell: Warning - ignoring invalid '
435 'unit specifier in %BLOCK {0} '
436 '(assuming Angstrom instead)'.format(blockname))
437 line_tokens = line_tokens[1:]
438 return u, line_tokens
440 # Arguments to pass to the Atoms object at the end
441 aargs = {
442 'pbc': True
443 }
445 # Start by looking for the lattice
446 lat_keywords = [w in celldict for w in ('lattice_cart', 'lattice_abc')]
447 if all(lat_keywords):
448 warnings.warn('read_cell: Warning - two lattice blocks present in the'
449 ' same file. LATTICE_ABC will be ignored')
450 elif not any(lat_keywords):
451 raise ValueError('Cell file must contain at least one between '
452 'LATTICE_ABC and LATTICE_CART')
454 if 'lattice_abc' in celldict:
456 lines = celldict.pop('lattice_abc')[0].split('\n')
457 line_tokens = [l.split() for l in lines]
459 u, line_tokens = parse_blockunit(line_tokens, 'lattice_abc')
461 if len(line_tokens) != 2:
462 warnings.warn('read_cell: Warning - ignoring additional '
463 'lines in invalid %BLOCK LATTICE_ABC')
465 abc = [float(p) * u for p in line_tokens[0][:3]]
466 angles = [float(phi) for phi in line_tokens[1][:3]]
468 aargs['cell'] = cellpar_to_cell(abc + angles)
470 if 'lattice_cart' in celldict:
472 lines = celldict.pop('lattice_cart')[0].split('\n')
473 line_tokens = [l.split() for l in lines]
475 u, line_tokens = parse_blockunit(line_tokens, 'lattice_cart')
477 if len(line_tokens) != 3:
478 warnings.warn('read_cell: Warning - ignoring more than '
479 'three lattice vectors in invalid %BLOCK '
480 'LATTICE_CART')
482 aargs['cell'] = [[float(x) * u for x in lt[:3]] for lt in line_tokens]
484 # Now move on to the positions
485 pos_keywords = [w in celldict
486 for w in ('positions_abs', 'positions_frac')]
488 if all(pos_keywords):
489 warnings.warn('read_cell: Warning - two lattice blocks present in the'
490 ' same file. POSITIONS_FRAC will be ignored')
491 del celldict['positions_frac']
492 elif not any(pos_keywords):
493 raise ValueError('Cell file must contain at least one between '
494 'POSITIONS_FRAC and POSITIONS_ABS')
496 aargs['symbols'] = []
497 pos_type = 'positions'
498 pos_block = celldict.pop('positions_abs', [None])[0]
499 if pos_block is None:
500 pos_type = 'scaled_positions'
501 pos_block = celldict.pop('positions_frac', [None])[0]
502 aargs[pos_type] = []
504 lines = pos_block.split('\n')
505 line_tokens = [l.split() for l in lines]
507 if 'scaled' not in pos_type:
508 u, line_tokens = parse_blockunit(line_tokens, 'positions_abs')
509 else:
510 u = 1.0
512 # Here we extract all the possible additional info
513 # These are marked by their type
515 add_info = {
516 'SPIN': (float, 0.0), # (type, default)
517 'MAGMOM': (float, 0.0),
518 'LABEL': (str, 'NULL')
519 }
520 add_info_arrays = dict((k, []) for k in add_info)
522 def parse_info(raw_info):
524 re_keys = (r'({0})\s*[=:\s]{{1}}\s'
525 r'*([^\s]*)').format('|'.join(add_info.keys()))
526 # Capture all info groups
527 info = re.findall(re_keys, raw_info)
528 info = {g[0]: add_info[g[0]][0](g[1]) for g in info}
529 return info
531 # Array for custom species (a CASTEP special thing)
532 # Usually left unused
533 custom_species = None
535 for tokens in line_tokens:
536 # Now, process the whole 'species' thing
537 spec_custom = tokens[0].split(':', 1)
538 elem = spec_custom[0]
539 if len(spec_custom) > 1 and custom_species is None:
540 # Add it to the custom info!
541 custom_species = list(aargs['symbols'])
542 if custom_species is not None:
543 custom_species.append(tokens[0])
544 aargs['symbols'].append(elem)
545 aargs[pos_type].append([float(p) * u for p in tokens[1:4]])
546 # Now for the additional information
547 info = ' '.join(tokens[4:])
548 info = parse_info(info)
549 for k in add_info:
550 add_info_arrays[k] += [info.get(k, add_info[k][1])]
552 # Now on to the species potentials...
553 if 'species_pot' in celldict:
554 lines = celldict.pop('species_pot')[0].split('\n')
555 line_tokens = [l.split() for l in lines]
557 for tokens in line_tokens:
558 if len(tokens) == 1:
559 # It's a library
560 all_spec = (set(custom_species) if custom_species is not None
561 else set(aargs['symbols']))
562 for s in all_spec:
563 calc.cell.species_pot = (s, tokens[0])
564 else:
565 calc.cell.species_pot = tuple(tokens[:2])
567 # Ionic constraints
568 raw_constraints = {}
570 if 'ionic_constraints' in celldict:
571 lines = celldict.pop('ionic_constraints')[0].split('\n')
572 line_tokens = [l.split() for l in lines]
574 for tokens in line_tokens:
575 if not len(tokens) == 6:
576 continue
577 _, species, nic, x, y, z = tokens
578 # convert xyz to floats
579 x = float(x)
580 y = float(y)
581 z = float(z)
583 nic = int(nic)
584 if (species, nic) not in raw_constraints:
585 raw_constraints[(species, nic)] = []
586 raw_constraints[(species, nic)].append(np.array(
587 [x, y, z]))
589 # Symmetry operations
590 if 'symmetry_ops' in celldict:
591 lines = celldict.pop('symmetry_ops')[0].split('\n')
592 line_tokens = [l.split() for l in lines]
594 # Read them in blocks of four
595 blocks = np.array(line_tokens).astype(float)
596 if (len(blocks.shape) != 2 or blocks.shape[1] != 3
597 or blocks.shape[0] % 4 != 0):
598 warnings.warn('Warning: could not parse SYMMETRY_OPS'
599 ' block properly, skipping')
600 else:
601 blocks = blocks.reshape((-1, 4, 3))
602 rotations = blocks[:, :3]
603 translations = blocks[:, 3]
605 # Regardless of whether we recognize them, store these
606 calc.cell.symmetry_ops = (rotations, translations)
608 # Anything else that remains, just add it to the cell object:
609 for k, (val, otype) in celldict.items():
610 try:
611 if otype == 'block':
612 val = val.split('\n') # Avoids a bug for one-line blocks
613 calc.cell.__setattr__(k, val)
614 except Exception as e:
615 raise RuntimeError('Problem setting calc.cell.%s = %s: %s' % (k, val, e))
617 # Get the relevant additional info
618 aargs['magmoms'] = np.array(add_info_arrays['SPIN'])
619 # SPIN or MAGMOM are alternative keywords
620 aargs['magmoms'] = np.where(aargs['magmoms'] != 0,
621 aargs['magmoms'],
622 add_info_arrays['MAGMOM'])
623 labels = np.array(add_info_arrays['LABEL'])
625 aargs['calculator'] = calc
627 atoms = ase.Atoms(**aargs)
629 # Spacegroup...
630 if find_spg:
631 # Try importing spglib
632 try:
633 import spglib
634 except ImportError:
635 warnings.warn('spglib not found installed on this system - '
636 'automatic spacegroup detection is not possible')
637 spglib = None
639 if spglib is not None:
640 symmd = spglib.get_symmetry_dataset(atoms_to_spglib_cell(atoms))
641 atoms_spg = Spacegroup(int(symmd['number']))
642 atoms.info['spacegroup'] = atoms_spg
644 atoms.new_array('castep_labels', labels)
645 if custom_species is not None:
646 atoms.new_array('castep_custom_species', np.array(custom_species))
648 fixed_atoms = []
649 constraints = []
650 for (species, nic), value in raw_constraints.items():
651 absolute_nr = atoms.calc._get_absolute_number(species, nic)
652 if len(value) == 3:
653 # Check if they are linearly independent
654 if np.linalg.det(value) == 0:
655 warnings.warn('Error: Found linearly dependent constraints attached '
656 'to atoms %s' % (absolute_nr))
657 continue
658 fixed_atoms.append(absolute_nr)
659 elif len(value) == 2:
660 direction = np.cross(value[0], value[1])
661 # Check if they are linearly independent
662 if np.linalg.norm(direction) == 0:
663 warnings.warn('Error: Found linearly dependent constraints attached '
664 'to atoms %s' % (absolute_nr))
665 continue
666 constraint = ase.constraints.FixedLine(
667 a=absolute_nr,
668 direction=direction)
669 constraints.append(constraint)
670 elif len(value) == 1:
671 constraint = ase.constraints.FixedPlane(
672 a=absolute_nr,
673 direction=np.array(value[0], dtype=np.float32))
674 constraints.append(constraint)
675 else:
676 warnings.warn('Error: Found %s statements attached to atoms %s' %
677 (len(value), absolute_nr))
679 # we need to sort the fixed atoms list in order not to raise an assertion
680 # error in FixAtoms
681 if fixed_atoms:
682 constraints.append(
683 ase.constraints.FixAtoms(indices=sorted(fixed_atoms)))
684 if constraints:
685 atoms.set_constraint(constraints)
687 atoms.calc.atoms = atoms
688 atoms.calc.push_oldstate()
690 return atoms
693def read_castep(filename, index=None):
694 """
695 Wrapper function for the more generic read() functionality.
697 Note that this is function is intended to maintain backwards-compatibility
698 only.
699 """
700 from ase.io import read
701 return read(filename, index=index, format='castep-castep')
704def read_castep_castep(fd, index=None):
705 """
706 Reads a .castep file and returns an atoms object.
707 The calculator information will be stored in the calc attribute.
709 There is no use of the "index" argument as of now, it is just inserted for
710 convenience to comply with the generic "read()" in ase.io
712 Please note that this routine will return an atom ordering as found
713 within the castep file. This means that the species will be ordered by
714 ascending atomic numbers. The atoms witin a species are ordered as given
715 in the original cell file.
717 Note: This routine returns a single atoms_object only, the last
718 configuration in the file. Yet, if you want to parse an MD run, use the
719 novel function `read_md()`
720 """
722 from ase.calculators.castep import Castep
724 try:
725 calc = Castep()
726 except Exception as e:
727 # No CASTEP keywords found?
728 warnings.warn('WARNING: {0} Using fallback .castep reader...'.format(e))
729 # Fall back on the old method
730 return read_castep_castep_old(fd, index)
732 calc.read(castep_file=fd)
734 # now we trick the calculator instance such that we can savely extract
735 # energies and forces from this atom. Basically what we do is to trick the
736 # internal routine calculation_required() to always return False such that
737 # we do not need to re-run a CASTEP calculation.
738 #
739 # Probably we can solve this with a flag to the read() routine at some
740 # point, but for the moment I do not want to change too much in there.
741 calc._old_atoms = calc.atoms
742 calc._old_param = calc.param
743 calc._old_cell = calc.cell
745 return [calc.atoms] # Returning in the form of a list for next()
748def read_castep_castep_old(fd, index=None):
749 """
750 DEPRECATED
751 Now replaced by ase.calculators.castep.Castep.read(). Left in for future
752 reference and backwards compatibility needs, as well as a fallback for
753 when castep_keywords.py can't be created.
755 Reads a .castep file and returns an atoms object.
756 The calculator information will be stored in the calc attribute.
757 If more than one SCF step is found, a list of all steps
758 will be stored in the traj attribute.
760 Note that the index argument has no effect as of now.
762 Please note that this routine will return an atom ordering as found
763 within the castep file. This means that the species will be ordered by
764 ascending atomic numbers. The atoms witin a species are ordered as given
765 in the original cell file.
766 """
767 from ase.calculators.singlepoint import SinglePointCalculator
769 lines = fd.readlines()
771 traj = []
772 energy_total = None
773 energy_0K = None
774 for i, line in enumerate(lines):
775 if 'NB est. 0K energy' in line:
776 energy_0K = float(line.split()[6])
777 # support also for dispersion correction
778 elif 'NB dispersion corrected est. 0K energy*' in line:
779 energy_0K = float(line.split()[-2])
780 elif 'Final energy, E' in line:
781 energy_total = float(line.split()[4])
782 elif 'Dispersion corrected final energy' in line:
783 pass
784 # dispcorr_energy_total = float(line.split()[-2])
785 # sedc_apply = True
786 elif 'Dispersion corrected final free energy' in line:
787 pass # dispcorr_energy_free = float(line.split()[-2])
788 elif 'dispersion corrected est. 0K energy' in line:
789 pass # dispcorr_energy_0K = float(line.split()[-2])
790 elif 'Unit Cell' in line:
791 cell = [x.split()[0:3] for x in lines[i + 3:i + 6]]
792 cell = np.array([[float(col) for col in row] for row in cell])
793 elif 'Cell Contents' in line:
794 geom_starts = i
795 start_found = False
796 for j, jline in enumerate(lines[geom_starts:]):
797 if jline.find('xxxxx') > 0 and start_found:
798 geom_stop = j + geom_starts
799 break
800 if jline.find('xxxx') > 0 and not start_found:
801 geom_start = j + geom_starts + 4
802 start_found = True
803 species = [line.split()[1] for line in lines[geom_start:geom_stop]]
804 geom = np.dot(np.array([[float(col) for col in line.split()[3:6]]
805 for line in lines[geom_start:geom_stop]]),
806 cell)
807 elif 'Writing model to' in line:
808 atoms = ase.Atoms(
809 cell=cell,
810 pbc=True,
811 positions=geom,
812 symbols=''.join(species))
813 # take 0K energy where available, else total energy
814 if energy_0K:
815 energy = energy_0K
816 else:
817 energy = energy_total
818 # generate a minimal single-point calculator
819 sp_calc = SinglePointCalculator(atoms=atoms,
820 energy=energy,
821 forces=None,
822 magmoms=None,
823 stress=None)
824 atoms.calc = sp_calc
825 traj.append(atoms)
826 if index is None:
827 return traj
828 else:
829 return traj[index]
832def read_geom(filename, index=':', units=units_CODATA2002):
833 """
834 Wrapper function for the more generic read() functionality.
836 Note that this is function is intended to maintain backwards-compatibility
837 only. Keyword arguments will be passed to read_castep_geom().
838 """
839 from ase.io import read
840 return read(filename, index=index, format='castep-geom', units=units)
843def read_castep_geom(fd, index=None, units=units_CODATA2002):
844 """Reads a .geom file produced by the CASTEP GeometryOptimization task and
845 returns an atoms object.
846 The information about total free energy and forces of each atom for every
847 relaxation step will be stored for further analysis especially in a
848 single-point calculator.
849 Note that everything in the .geom file is in atomic units, which has
850 been conversed to commonly used unit angstrom(length) and eV (energy).
852 Note that the index argument has no effect as of now.
854 Contribution by Wei-Bing Zhang. Thanks!
856 Routine now accepts a filedescriptor in order to out-source the gz and
857 bz2 handling to formats.py. Note that there is a fallback routine
858 read_geom() that behaves like previous versions did.
859 """
860 from ase.calculators.singlepoint import SinglePointCalculator
862 # fd is closed by embracing read() routine
863 txt = fd.readlines()
865 traj = []
867 Hartree = units['Eh']
868 Bohr = units['a0']
870 # Yeah, we know that...
871 # print('N.B.: Energy in .geom file is not 0K extrapolated.')
872 for i, line in enumerate(txt):
873 if line.find('<-- E') > 0:
874 start_found = True
875 energy = float(line.split()[0]) * Hartree
876 cell = [x.split()[0:3] for x in txt[i + 1:i + 4]]
877 cell = np.array([[float(col) * Bohr for col in row] for row in
878 cell])
879 if line.find('<-- R') > 0 and start_found:
880 start_found = False
881 geom_start = i
882 for i, line in enumerate(txt[geom_start:]):
883 if line.find('<-- F') > 0:
884 geom_stop = i + geom_start
885 break
886 species = [line.split()[0] for line in
887 txt[geom_start:geom_stop]]
888 geom = np.array([[float(col) * Bohr for col in
889 line.split()[2:5]] for line in
890 txt[geom_start:geom_stop]])
891 forces = np.array([[float(col) * Hartree / Bohr for col in
892 line.split()[2:5]] for line in
893 txt[geom_stop:geom_stop
894 + (geom_stop - geom_start)]])
895 image = ase.Atoms(species, geom, cell=cell, pbc=True)
896 image.calc = SinglePointCalculator(
897 atoms=image, energy=energy, forces=forces)
898 traj.append(image)
900 if index is None:
901 return traj
902 else:
903 return traj[index]
906def read_phonon(filename, index=None, read_vib_data=False,
907 gamma_only=True, frequency_factor=None,
908 units=units_CODATA2002):
909 """
910 Wrapper function for the more generic read() functionality.
912 Note that this is function is intended to maintain backwards-compatibility
913 only. For documentation see read_castep_phonon().
914 """
915 from ase.io import read
917 if read_vib_data:
918 full_output = True
919 else:
920 full_output = False
922 return read(filename, index=index, format='castep-phonon',
923 full_output=full_output, read_vib_data=read_vib_data,
924 gamma_only=gamma_only, frequency_factor=frequency_factor,
925 units=units)
928def read_castep_phonon(fd, index=None, read_vib_data=False,
929 gamma_only=True, frequency_factor=None,
930 units=units_CODATA2002):
931 """
932 Reads a .phonon file written by a CASTEP Phonon task and returns an atoms
933 object, as well as the calculated vibrational data if requested.
935 Note that the index argument has no effect as of now.
936 """
938 # fd is closed by embracing read() routine
939 lines = fd.readlines()
941 atoms = None
942 cell = []
943 N = Nb = Nq = 0
944 scaled_positions = []
945 symbols = []
946 masses = []
948 # header
949 L = 0
950 while L < len(lines):
952 line = lines[L]
954 if 'Number of ions' in line:
955 N = int(line.split()[3])
956 elif 'Number of branches' in line:
957 Nb = int(line.split()[3])
958 elif 'Number of wavevectors' in line:
959 Nq = int(line.split()[3])
960 elif 'Unit cell vectors (A)' in line:
961 for ll in range(3):
962 L += 1
963 fields = lines[L].split()
964 cell.append([float(x) for x in fields[0:3]])
965 elif 'Fractional Co-ordinates' in line:
966 for ll in range(N):
967 L += 1
968 fields = lines[L].split()
969 scaled_positions.append([float(x) for x in fields[1:4]])
970 symbols.append(fields[4])
971 masses.append(float(fields[5]))
972 elif 'END header' in line:
973 L += 1
974 atoms = ase.Atoms(symbols=symbols,
975 scaled_positions=scaled_positions,
976 cell=cell)
977 break
979 L += 1
981 # Eigenmodes and -vectors
982 if frequency_factor is None:
983 Kayser_to_eV = 1E2 * 2 * np.pi * units['hbar'] * units['c']
984 # N.B. "fixed default" unit for frequencies in .phonon files is "cm-1"
985 # (i.e. the latter is unaffected by the internal unit conversion system of
986 # CASTEP!) set conversion factor to convert therefrom to eV by default for
987 # now
988 frequency_factor = Kayser_to_eV
989 qpoints = []
990 weights = []
991 frequencies = []
992 displacements = []
993 for nq in range(Nq):
994 fields = lines[L].split()
995 qpoints.append([float(x) for x in fields[2:5]])
996 weights.append(float(fields[5]))
997 freqs = []
998 for ll in range(Nb):
999 L += 1
1000 fields = lines[L].split()
1001 freqs.append(frequency_factor * float(fields[1]))
1002 frequencies.append(np.array(freqs))
1004 # skip the two Phonon Eigenvectors header lines
1005 L += 2
1007 # generate a list of displacements with a structure that is identical to
1008 # what is stored internally in the Vibrations class (see in
1009 # ase.vibrations.Vibrations.modes):
1010 # np.array(displacements).shape == (Nb,3*N)
1012 disps = []
1013 for ll in range(Nb):
1014 disp_coords = []
1015 for lll in range(N):
1016 L += 1
1017 fields = lines[L].split()
1018 disp_x = float(fields[2]) + float(fields[3]) * 1.0j
1019 disp_y = float(fields[4]) + float(fields[5]) * 1.0j
1020 disp_z = float(fields[6]) + float(fields[7]) * 1.0j
1021 disp_coords.extend([disp_x, disp_y, disp_z])
1022 disps.append(np.array(disp_coords))
1023 displacements.append(np.array(disps))
1025 if read_vib_data:
1026 if gamma_only:
1027 vibdata = [frequencies[0], displacements[0]]
1028 else:
1029 vibdata = [qpoints, weights, frequencies, displacements]
1030 return vibdata, atoms
1031 else:
1032 return atoms
1035def read_md(filename, index=None, return_scalars=False,
1036 units=units_CODATA2002):
1037 """Wrapper function for the more generic read() functionality.
1039 Note that this function is intended to maintain backwards-compatibility
1040 only. For documentation see read_castep_md()
1041 """
1042 if return_scalars:
1043 full_output = True
1044 else:
1045 full_output = False
1047 from ase.io import read
1048 return read(filename, index=index, format='castep-md',
1049 full_output=full_output, return_scalars=return_scalars,
1050 units=units)
1053def read_castep_md(fd, index=None, return_scalars=False,
1054 units=units_CODATA2002):
1055 """Reads a .md file written by a CASTEP MolecularDynamics task
1056 and returns the trajectory stored therein as a list of atoms object.
1058 Note that the index argument has no effect as of now."""
1060 from ase.calculators.singlepoint import SinglePointCalculator
1062 factors = {
1063 't': units['t0'] * 1E15, # fs
1064 'E': units['Eh'], # eV
1065 'T': units['Eh'] / units['kB'],
1066 'P': units['Eh'] / units['a0']**3 * units['Pascal'],
1067 'h': units['a0'],
1068 'hv': units['a0'] / units['t0'],
1069 'S': units['Eh'] / units['a0']**3,
1070 'R': units['a0'],
1071 'V': np.sqrt(units['Eh'] / units['me']),
1072 'F': units['Eh'] / units['a0']}
1074 # fd is closed by embracing read() routine
1075 lines = fd.readlines()
1077 L = 0
1078 while 'END header' not in lines[L]:
1079 L += 1
1080 l_end_header = L
1081 lines = lines[l_end_header + 1:]
1082 times = []
1083 energies = []
1084 temperatures = []
1085 pressures = []
1086 traj = []
1088 # Initialization
1089 time = None
1090 Epot = None
1091 Ekin = None
1092 EH = None
1093 temperature = None
1094 pressure = None
1095 symbols = None
1096 positions = None
1097 cell = None
1098 velocities = None
1099 symbols = []
1100 positions = []
1101 velocities = []
1102 forces = []
1103 cell = np.eye(3)
1104 cell_velocities = []
1105 stress = []
1107 for (L, line) in enumerate(lines):
1108 fields = line.split()
1109 if len(fields) == 0:
1110 if L != 0:
1111 times.append(time)
1112 energies.append([Epot, EH, Ekin])
1113 temperatures.append(temperature)
1114 pressures.append(pressure)
1115 atoms = ase.Atoms(symbols=symbols,
1116 positions=positions,
1117 cell=cell)
1118 atoms.set_velocities(velocities)
1119 if len(stress) == 0:
1120 atoms.calc = SinglePointCalculator(
1121 atoms=atoms, energy=Epot, forces=forces)
1122 else:
1123 atoms.calc = SinglePointCalculator(
1124 atoms=atoms, energy=Epot,
1125 forces=forces, stress=stress)
1126 traj.append(atoms)
1127 symbols = []
1128 positions = []
1129 velocities = []
1130 forces = []
1131 cell = []
1132 cell_velocities = []
1133 stress = []
1134 continue
1135 if len(fields) == 1:
1136 time = factors['t'] * float(fields[0])
1137 continue
1139 if fields[-1] == 'E':
1140 E = [float(x) for x in fields[0:3]]
1141 Epot, EH, Ekin = [factors['E'] * Ei for Ei in E]
1142 continue
1144 if fields[-1] == 'T':
1145 temperature = factors['T'] * float(fields[0])
1146 continue
1148 # only printed in case of variable cell calculation or calculate_stress
1149 # explicitly requested
1150 if fields[-1] == 'P':
1151 pressure = factors['P'] * float(fields[0])
1152 continue
1153 if fields[-1] == 'h':
1154 h = [float(x) for x in fields[0:3]]
1155 cell.append([factors['h'] * hi for hi in h])
1156 continue
1158 # only printed in case of variable cell calculation
1159 if fields[-1] == 'hv':
1160 hv = [float(x) for x in fields[0:3]]
1161 cell_velocities.append([factors['hv'] * hvi for hvi in hv])
1162 continue
1164 # only printed in case of variable cell calculation
1165 if fields[-1] == 'S':
1166 S = [float(x) for x in fields[0:3]]
1167 stress.append([factors['S'] * Si for Si in S])
1168 continue
1169 if fields[-1] == 'R':
1170 symbols.append(fields[0])
1171 R = [float(x) for x in fields[2:5]]
1172 positions.append([factors['R'] * Ri for Ri in R])
1173 continue
1174 if fields[-1] == 'V':
1175 V = [float(x) for x in fields[2:5]]
1176 velocities.append([factors['V'] * Vi for Vi in V])
1177 continue
1178 if fields[-1] == 'F':
1179 F = [float(x) for x in fields[2:5]]
1180 forces.append([factors['F'] * Fi for Fi in F])
1181 continue
1183 if index is None:
1184 pass
1185 else:
1186 traj = traj[index]
1188 if return_scalars:
1189 data = [times, energies, temperatures, pressures]
1190 return data, traj
1191 else:
1192 return traj
1195# Routines that only the calculator requires
1197def read_param(filename='', calc=None, fd=None, get_interface_options=False):
1199 if fd is None:
1200 if filename == '':
1201 raise ValueError('One between filename and fd must be provided')
1202 fd = open(filename)
1203 elif filename:
1204 warnings.warn('Filestream used to read param, file name will be '
1205 'ignored')
1207 # If necessary, get the interface options
1208 if get_interface_options:
1209 int_opts = {}
1210 optre = re.compile(r'# ASE_INTERFACE ([^\s]+) : ([^\s]+)')
1212 lines = fd.readlines()
1213 fd.seek(0)
1215 for l in lines:
1216 m = optre.search(l)
1217 if m:
1218 int_opts[m.groups()[0]] = m.groups()[1]
1220 data = read_freeform(fd)
1222 if calc is None:
1223 from ase.calculators.castep import Castep
1224 calc = Castep(check_castep_version=False, keyword_tolerance=2)
1226 for kw, (val, otype) in data.items():
1227 if otype == 'block':
1228 val = val.split('\n') # Avoids a bug for one-line blocks
1229 calc.param.__setattr__(kw, val)
1231 if not get_interface_options:
1232 return calc
1233 else:
1234 return calc, int_opts
1237def write_param(filename, param, check_checkfile=False,
1238 force_write=False,
1239 interface_options=None):
1240 """Writes a CastepParam object to a CASTEP .param file
1242 Parameters:
1243 filename: the location of the file to write to. If it
1244 exists it will be overwritten without warning. If it
1245 doesn't it will be created.
1246 param: a CastepParam instance
1247 check_checkfile : if set to True, write_param will
1248 only write continuation or reuse statement
1249 if a restart file exists in the same directory
1250 """
1251 if os.path.isfile(filename) and not force_write:
1252 warnings.warn('ase.io.castep.write_param: Set optional argument '
1253 'force_write=True to overwrite %s.' % filename)
1254 return False
1256 out = paropen(filename, 'w')
1257 out.write('#######################################################\n')
1258 out.write('#CASTEP param file: %s\n' % filename)
1259 out.write('#Created using the Atomic Simulation Environment (ASE)#\n')
1260 if interface_options is not None:
1261 out.write('# Internal settings of the calculator\n')
1262 out.write('# This can be switched off by settings\n')
1263 out.write('# calc._export_settings = False\n')
1264 out.write('# If stated, this will be automatically processed\n')
1265 out.write('# by ase.io.castep.read_seed()\n')
1266 for option, value in sorted(interface_options.items()):
1267 out.write('# ASE_INTERFACE %s : %s\n' % (option, value))
1268 out.write('#######################################################\n\n')
1270 if check_checkfile:
1271 param = deepcopy(param) # To avoid modifying the parent one
1272 for checktype in ['continuation', 'reuse']:
1273 opt = getattr(param, checktype)
1274 if opt and opt.value:
1275 fname = opt.value
1276 if fname == 'default':
1277 fname = os.path.splitext(filename)[0] + '.check'
1278 if not (os.path.exists(fname) or
1279 # CASTEP also understands relative path names, hence
1280 # also check relative to the param file directory
1281 os.path.exists(
1282 os.path.join(os.path.dirname(filename),
1283 opt.value))):
1284 opt.clear()
1286 write_freeform(out, param)
1288 out.close()
1291def read_seed(seed, new_seed=None, ignore_internal_keys=False):
1292 """A wrapper around the CASTEP Calculator in conjunction with
1293 read_cell and read_param. Basically this can be used to reuse
1294 a previous calculation which results in a triple of
1295 cell/param/castep file. The label of the calculation if pre-
1296 fixed with `copy_of_` and everything else will be recycled as
1297 much as possible from the addressed calculation.
1299 Please note that this routine will return an atoms ordering as specified
1300 in the cell file! It will thus undo the potential reordering internally
1301 done by castep.
1302 """
1304 directory = os.path.abspath(os.path.dirname(seed))
1305 seed = os.path.basename(seed)
1307 paramfile = os.path.join(directory, '%s.param' % seed)
1308 cellfile = os.path.join(directory, '%s.cell' % seed)
1309 castepfile = os.path.join(directory, '%s.castep' % seed)
1310 checkfile = os.path.join(directory, '%s.check' % seed)
1312 atoms = read_cell(cellfile)
1313 atoms.calc._directory = directory
1314 atoms.calc._rename_existing_dir = False
1315 atoms.calc._castep_pp_path = directory
1316 atoms.calc.merge_param(paramfile,
1317 ignore_internal_keys=ignore_internal_keys)
1318 if new_seed is None:
1319 atoms.calc._label = 'copy_of_%s' % seed
1320 else:
1321 atoms.calc._label = str(new_seed)
1322 if os.path.isfile(castepfile):
1323 # _set_atoms needs to be True here
1324 # but we set it right back to False
1325 # atoms.calc._set_atoms = False
1326 # BUGFIX: I do not see a reason to do that!
1327 atoms.calc.read(castepfile)
1328 # atoms.calc._set_atoms = False
1330 # if here is a check file, we also want to re-use this information
1331 if os.path.isfile(checkfile):
1332 atoms.calc._check_file = os.path.basename(checkfile)
1334 # sync the top-level object with the
1335 # one attached to the calculator
1336 atoms = atoms.calc.atoms
1337 else:
1338 # There are cases where we only want to restore a calculator/atoms
1339 # setting without a castep file...
1340 pass
1341 # No print statement required in these cases
1342 warnings.warn('Corresponding *.castep file not found. '
1343 'Atoms object will be restored from *.cell and *.param only.')
1344 atoms.calc.push_oldstate()
1346 return atoms
1349def read_bands(filename='', fd=None, units=units_CODATA2002):
1350 """Read Castep.bands file to kpoints, weights and eigenvalues
1352 Args:
1353 filename (str):
1354 path to seedname.bands file
1355 fd (fd):
1356 file descriptor for open bands file
1357 units (dict):
1358 Conversion factors for atomic units
1360 Returns:
1361 (tuple):
1362 (kpts, weights, eigenvalues, efermi)
1364 Where ``kpts`` and ``weights`` are 1d numpy arrays, eigenvalues
1365 is an array of shape (spin, kpts, nbands) and efermi is a float
1366 """
1368 Hartree = units['Eh']
1370 if fd is None:
1371 if filename == '':
1372 raise ValueError('One between filename and fd must be provided')
1373 fd = open(filename)
1374 elif filename:
1375 warnings.warn('Filestream used to read param, file name will be '
1376 'ignored')
1378 nkpts, nspin, _, nbands, efermi = [t(fd.readline().split()[-1]) for t in
1379 [int, int, float, int, float]]
1381 kpts, weights = np.zeros((nkpts, 3)), np.zeros(nkpts)
1382 eigenvalues = np.zeros((nspin, nkpts, nbands))
1384 # Skip unit cell
1385 for _ in range(4):
1386 fd.readline()
1388 def _kptline_to_i_k_wt(line):
1389 line = line.split()
1390 line = [int(line[1])] + list(map(float, line[2:]))
1391 return (line[0] - 1, line[1:4], line[4])
1393 # CASTEP often writes these out-of-order, so check index and write directly
1394 # to the correct row
1395 for kpt_line in range(nkpts):
1396 i_kpt, kpt, wt = _kptline_to_i_k_wt(fd.readline())
1397 kpts[i_kpt, :], weights[i_kpt] = kpt, wt
1398 for spin in range(nspin):
1399 fd.readline() # Skip 'Spin component N' line
1400 eigenvalues[spin, i_kpt, :] = [float(fd.readline())
1401 for _ in range(nbands)]
1403 return (kpts, weights, eigenvalues * Hartree, efermi * Hartree)