Coverage for /builds/debichem-team/python-ase/ase/io/vasp.py : 77.83%
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"""
2This module contains functionality for reading and writing an ASE
3Atoms object in VASP POSCAR format.
5"""
7import re
9import numpy as np
11from ase import Atoms
12from ase.utils import reader, writer
13from ase.io.utils import ImageIterator
14from ase.io import ParseError
15from .vasp_parsers import vasp_outcar_parsers as vop
16from pathlib import Path
18__all__ = [
19 'read_vasp', 'read_vasp_out', 'iread_vasp_out', 'read_vasp_xdatcar',
20 'read_vasp_xml', 'write_vasp', 'write_vasp_xdatcar'
21]
24def get_atomtypes(fname):
25 """Given a file name, get the atomic symbols.
27 The function can get this information from OUTCAR and POTCAR
28 format files. The files can also be compressed with gzip or
29 bzip2.
31 """
32 fpath = Path(fname)
34 atomtypes = []
35 atomtypes_alt = []
36 if fpath.suffix == '.gz':
37 import gzip
38 opener = gzip.open
39 elif fpath.suffix == '.bz2':
40 import bz2
41 opener = bz2.BZ2File
42 else:
43 opener = open
44 with opener(fpath) as fd:
45 for line in fd:
46 if 'TITEL' in line:
47 atomtypes.append(line.split()[3].split('_')[0].split('.')[0])
48 elif 'POTCAR:' in line:
49 atomtypes_alt.append(line.split()[2].split('_')[0].split('.')[0])
51 if len(atomtypes) == 0 and len(atomtypes_alt) > 0:
52 # old VASP doesn't echo TITEL, but all versions print out species lines
53 # preceded by "POTCAR:", twice
54 if len(atomtypes_alt) % 2 != 0:
55 raise ParseError(f'Tried to get atom types from {len(atomtypes_alt)} "POTCAR": '
56 'lines in OUTCAR, but expected an even number')
57 atomtypes = atomtypes_alt[0:len(atomtypes_alt)//2]
59 return atomtypes
62def atomtypes_outpot(posfname, numsyms):
63 """Try to retrieve chemical symbols from OUTCAR or POTCAR
65 If getting atomtypes from the first line in POSCAR/CONTCAR fails, it might
66 be possible to find the data in OUTCAR or POTCAR, if these files exist.
68 posfname -- The filename of the POSCAR/CONTCAR file we're trying to read
70 numsyms -- The number of symbols we must find
72 """
73 posfpath = Path(posfname)
75 # Check files with exactly same path except POTCAR/OUTCAR instead
76 # of POSCAR/CONTCAR.
77 fnames = [posfpath.with_name('POTCAR'),
78 posfpath.with_name('OUTCAR')]
79 # Try the same but with compressed files
80 fsc = []
81 for fnpath in fnames:
82 fsc.append(fnpath.parent / (fnpath.name + '.gz'))
83 fsc.append(fnpath.parent / (fnpath.name + '.bz2'))
84 for f in fsc:
85 fnames.append(f)
86 # Code used to try anything with POTCAR or OUTCAR in the name
87 # but this is no longer supported
89 tried = []
90 for fn in fnames:
91 if fn in posfpath.parent.iterdir():
92 tried.append(fn)
93 at = get_atomtypes(fn)
94 if len(at) == numsyms:
95 return at
97 raise ParseError('Could not determine chemical symbols. Tried files ' +
98 str(tried))
101def get_atomtypes_from_formula(formula):
102 """Return atom types from chemical formula (optionally prepended
103 with and underscore).
104 """
105 from ase.symbols import string2symbols
106 symbols = string2symbols(formula.split('_')[0])
107 atomtypes = [symbols[0]]
108 for s in symbols[1:]:
109 if s != atomtypes[-1]:
110 atomtypes.append(s)
111 return atomtypes
114@reader
115def read_vasp(filename='CONTCAR'):
116 """Import POSCAR/CONTCAR type file.
118 Reads unitcell, atom positions and constraints from the POSCAR/CONTCAR
119 file and tries to read atom types from POSCAR/CONTCAR header, if this fails
120 the atom types are read from OUTCAR or POTCAR file.
121 """
123 from ase.constraints import FixAtoms, FixScaled
124 from ase.data import chemical_symbols
126 fd = filename
127 # The first line is in principle a comment line, however in VASP
128 # 4.x a common convention is to have it contain the atom symbols,
129 # eg. "Ag Ge" in the same order as later in the file (and POTCAR
130 # for the full vasp run). In the VASP 5.x format this information
131 # is found on the fifth line. Thus we save the first line and use
132 # it in case we later detect that we're reading a VASP 4.x format
133 # file.
134 line1 = fd.readline()
136 lattice_constant = float(fd.readline().split()[0])
138 # Now the lattice vectors
139 a = []
140 for ii in range(3):
141 s = fd.readline().split()
142 floatvect = float(s[0]), float(s[1]), float(s[2])
143 a.append(floatvect)
145 basis_vectors = np.array(a) * lattice_constant
147 # Number of atoms. Again this must be in the same order as
148 # in the first line
149 # or in the POTCAR or OUTCAR file
150 atom_symbols = []
151 numofatoms = fd.readline().split()
152 # Check whether we have a VASP 4.x or 5.x format file. If the
153 # format is 5.x, use the fifth line to provide information about
154 # the atomic symbols.
155 vasp5 = False
156 try:
157 int(numofatoms[0])
158 except ValueError:
159 vasp5 = True
160 atomtypes = numofatoms
161 numofatoms = fd.readline().split()
163 # check for comments in numofatoms line and get rid of them if necessary
164 commentcheck = np.array(['!' in s for s in numofatoms])
165 if commentcheck.any():
166 # only keep the elements up to the first including a '!':
167 numofatoms = numofatoms[:np.arange(len(numofatoms))[commentcheck][0]]
169 if not vasp5:
170 # Split the comment line (first in the file) into words and
171 # try to compose a list of chemical symbols
172 from ase.formula import Formula
173 atomtypes = []
174 for word in line1.split():
175 word_without_delims = re.sub(r"-|_|,|\.|=|[0-9]|^", "", word)
176 if len(word_without_delims) < 1:
177 continue
178 try:
179 atomtypes.extend(list(Formula(word_without_delims)))
180 except ValueError:
181 # print(atomtype, e, 'is comment')
182 pass
183 # Now the list of chemical symbols atomtypes must be formed.
184 # For example: atomtypes = ['Pd', 'C', 'O']
186 numsyms = len(numofatoms)
187 if len(atomtypes) < numsyms:
188 # First line in POSCAR/CONTCAR didn't contain enough symbols.
190 # Sometimes the first line in POSCAR/CONTCAR is of the form
191 # "CoP3_In-3.pos". Check for this case and extract atom types
192 if len(atomtypes) == 1 and '_' in atomtypes[0]:
193 atomtypes = get_atomtypes_from_formula(atomtypes[0])
194 else:
195 atomtypes = atomtypes_outpot(fd.name, numsyms)
196 else:
197 try:
198 for atype in atomtypes[:numsyms]:
199 if atype not in chemical_symbols:
200 raise KeyError
201 except KeyError:
202 atomtypes = atomtypes_outpot(fd.name, numsyms)
204 for i, num in enumerate(numofatoms):
205 numofatoms[i] = int(num)
206 [atom_symbols.append(atomtypes[i]) for na in range(numofatoms[i])]
208 # Check if Selective dynamics is switched on
209 sdyn = fd.readline()
210 selective_dynamics = sdyn[0].lower() == 's'
212 # Check if atom coordinates are cartesian or direct
213 if selective_dynamics:
214 ac_type = fd.readline()
215 else:
216 ac_type = sdyn
217 cartesian = ac_type[0].lower() == 'c' or ac_type[0].lower() == 'k'
218 tot_natoms = sum(numofatoms)
219 atoms_pos = np.empty((tot_natoms, 3))
220 if selective_dynamics:
221 selective_flags = np.empty((tot_natoms, 3), dtype=bool)
222 for atom in range(tot_natoms):
223 ac = fd.readline().split()
224 atoms_pos[atom] = (float(ac[0]), float(ac[1]), float(ac[2]))
225 if selective_dynamics:
226 curflag = []
227 for flag in ac[3:6]:
228 curflag.append(flag == 'F')
229 selective_flags[atom] = curflag
230 if cartesian:
231 atoms_pos *= lattice_constant
232 atoms = Atoms(symbols=atom_symbols, cell=basis_vectors, pbc=True)
233 if cartesian:
234 atoms.set_positions(atoms_pos)
235 else:
236 atoms.set_scaled_positions(atoms_pos)
237 if selective_dynamics:
238 constraints = []
239 indices = []
240 for ind, sflags in enumerate(selective_flags):
241 if sflags.any() and not sflags.all():
242 constraints.append(FixScaled(atoms.get_cell(), ind, sflags))
243 elif sflags.all():
244 indices.append(ind)
245 if indices:
246 constraints.append(FixAtoms(indices))
247 if constraints:
248 atoms.set_constraint(constraints)
249 return atoms
252def iread_vasp_out(filename, index=-1):
253 """Import OUTCAR type file, as a generator."""
254 it = ImageIterator(vop.outcarchunks)
255 return it(filename, index=index)
258@reader
259def read_vasp_out(filename='OUTCAR', index=-1):
260 """Import OUTCAR type file.
262 Reads unitcell, atom positions, energies, and forces from the OUTCAR file
263 and attempts to read constraints (if any) from CONTCAR/POSCAR, if present.
264 """
265 # "filename" is actually a file-descriptor thanks to @reader
266 g = iread_vasp_out(filename, index=index)
267 # Code borrowed from formats.py:read
268 if isinstance(index, (slice, str)):
269 # Return list of atoms
270 return list(g)
271 else:
272 # Return single atoms object
273 return next(g)
276@reader
277def read_vasp_xdatcar(filename='XDATCAR', index=-1):
278 """Import XDATCAR file
280 Reads all positions from the XDATCAR and returns a list of
281 Atoms objects. Useful for viewing optimizations runs
282 from VASP5.x
284 Constraints ARE NOT stored in the XDATCAR, and as such, Atoms
285 objects retrieved from the XDATCAR will not have constraints set.
286 """
287 fd = filename # @reader decorator ensures this is a file descriptor
288 images = list()
290 cell = np.eye(3)
291 atomic_formula = str()
293 while True:
294 comment_line = fd.readline()
295 if "Direct configuration=" not in comment_line:
296 try:
297 lattice_constant = float(fd.readline())
298 except Exception:
299 # XXX: When would this happen?
300 break
302 xx = [float(x) for x in fd.readline().split()]
303 yy = [float(y) for y in fd.readline().split()]
304 zz = [float(z) for z in fd.readline().split()]
305 cell = np.array([xx, yy, zz]) * lattice_constant
307 symbols = fd.readline().split()
308 numbers = [int(n) for n in fd.readline().split()]
309 total = sum(numbers)
311 atomic_formula = ''.join('{:s}{:d}'.format(sym, numbers[n])
312 for n, sym in enumerate(symbols))
314 fd.readline()
316 coords = [
317 np.array(fd.readline().split(), float) for ii in range(total)
318 ]
320 image = Atoms(atomic_formula, cell=cell, pbc=True)
321 image.set_scaled_positions(np.array(coords))
322 images.append(image)
324 if not index:
325 return images
326 else:
327 return images[index]
330def __get_xml_parameter(par):
331 """An auxiliary function that enables convenient extraction of
332 parameter values from a vasprun.xml file with proper type
333 handling.
335 """
336 def to_bool(b):
337 if b == 'T':
338 return True
339 else:
340 return False
342 to_type = {'int': int, 'logical': to_bool, 'string': str, 'float': float}
344 text = par.text
345 if text is None:
346 text = ''
348 # Float parameters do not have a 'type' attrib
349 var_type = to_type[par.attrib.get('type', 'float')]
351 try:
352 if par.tag == 'v':
353 return list(map(var_type, text.split()))
354 else:
355 return var_type(text.strip())
356 except ValueError:
357 # Vasp can sometimes write "*****" due to overflow
358 return None
361def read_vasp_xml(filename='vasprun.xml', index=-1):
362 """Parse vasprun.xml file.
364 Reads unit cell, atom positions, energies, forces, and constraints
365 from vasprun.xml file
366 """
368 import xml.etree.ElementTree as ET
369 from ase.constraints import FixAtoms, FixScaled
370 from ase.calculators.singlepoint import (SinglePointDFTCalculator,
371 SinglePointKPoint)
372 from ase.units import GPa
373 from collections import OrderedDict
375 tree = ET.iterparse(filename, events=['start', 'end'])
377 atoms_init = None
378 calculation = []
379 ibz_kpts = None
380 kpt_weights = None
381 parameters = OrderedDict()
383 try:
384 for event, elem in tree:
386 if event == 'end':
387 if elem.tag == 'kpoints':
388 for subelem in elem.iter(tag='generation'):
389 kpts_params = OrderedDict()
390 parameters['kpoints_generation'] = kpts_params
391 for par in subelem.iter():
392 if par.tag in ['v', 'i']:
393 parname = par.attrib['name'].lower()
394 kpts_params[parname] = __get_xml_parameter(par)
396 kpts = elem.findall("varray[@name='kpointlist']/v")
397 ibz_kpts = np.zeros((len(kpts), 3))
399 for i, kpt in enumerate(kpts):
400 ibz_kpts[i] = [float(val) for val in kpt.text.split()]
402 kpt_weights = elem.findall('varray[@name="weights"]/v')
403 kpt_weights = [float(val.text) for val in kpt_weights]
405 elif elem.tag == 'parameters':
406 for par in elem.iter():
407 if par.tag in ['v', 'i']:
408 parname = par.attrib['name'].lower()
409 parameters[parname] = __get_xml_parameter(par)
411 elif elem.tag == 'atominfo':
412 species = []
414 for entry in elem.find("array[@name='atoms']/set"):
415 species.append(entry[0].text.strip())
417 natoms = len(species)
419 elif (elem.tag == 'structure'
420 and elem.attrib.get('name') == 'initialpos'):
421 cell_init = np.zeros((3, 3), dtype=float)
423 for i, v in enumerate(
424 elem.find("crystal/varray[@name='basis']")):
425 cell_init[i] = np.array(
426 [float(val) for val in v.text.split()])
428 scpos_init = np.zeros((natoms, 3), dtype=float)
430 for i, v in enumerate(
431 elem.find("varray[@name='positions']")):
432 scpos_init[i] = np.array(
433 [float(val) for val in v.text.split()])
435 constraints = []
436 fixed_indices = []
438 for i, entry in enumerate(
439 elem.findall("varray[@name='selective']/v")):
440 flags = (np.array(
441 entry.text.split() == np.array(['F', 'F', 'F'])))
442 if flags.all():
443 fixed_indices.append(i)
444 elif flags.any():
445 constraints.append(FixScaled(cell_init, i, flags))
447 if fixed_indices:
448 constraints.append(FixAtoms(fixed_indices))
450 atoms_init = Atoms(species,
451 cell=cell_init,
452 scaled_positions=scpos_init,
453 constraint=constraints,
454 pbc=True)
456 elif elem.tag == 'dipole':
457 dblock = elem.find('v[@name="dipole"]')
458 if dblock is not None:
459 dipole = np.array(
460 [float(val) for val in dblock.text.split()])
462 elif event == 'start' and elem.tag == 'calculation':
463 calculation.append(elem)
465 except ET.ParseError as parse_error:
466 if atoms_init is None:
467 raise parse_error
468 if calculation and calculation[-1].find("energy") is None:
469 calculation = calculation[:-1]
470 if not calculation:
471 yield atoms_init
473 if calculation:
474 if isinstance(index, int):
475 steps = [calculation[index]]
476 else:
477 steps = calculation[index]
478 else:
479 steps = []
481 for step in steps:
482 # Workaround for VASP bug, e_0_energy contains the wrong value
483 # in calculation/energy, but calculation/scstep/energy does not
484 # include classical VDW corrections. So, first calculate
485 # e_0_energy - e_fr_energy from calculation/scstep/energy, then
486 # apply that correction to e_fr_energy from calculation/energy.
487 lastscf = step.findall('scstep/energy')[-1]
488 dipoles = step.findall('scstep/dipole')
489 if dipoles:
490 lastdipole = dipoles[-1]
491 else:
492 lastdipole = None
494 de = (float(lastscf.find('i[@name="e_0_energy"]').text) -
495 float(lastscf.find('i[@name="e_fr_energy"]').text))
497 free_energy = float(step.find('energy/i[@name="e_fr_energy"]').text)
498 energy = free_energy + de
500 cell = np.zeros((3, 3), dtype=float)
501 for i, vector in enumerate(
502 step.find('structure/crystal/varray[@name="basis"]')):
503 cell[i] = np.array([float(val) for val in vector.text.split()])
505 scpos = np.zeros((natoms, 3), dtype=float)
506 for i, vector in enumerate(
507 step.find('structure/varray[@name="positions"]')):
508 scpos[i] = np.array([float(val) for val in vector.text.split()])
510 forces = None
511 fblocks = step.find('varray[@name="forces"]')
512 if fblocks is not None:
513 forces = np.zeros((natoms, 3), dtype=float)
514 for i, vector in enumerate(fblocks):
515 forces[i] = np.array(
516 [float(val) for val in vector.text.split()])
518 stress = None
519 sblocks = step.find('varray[@name="stress"]')
520 if sblocks is not None:
521 stress = np.zeros((3, 3), dtype=float)
522 for i, vector in enumerate(sblocks):
523 stress[i] = np.array(
524 [float(val) for val in vector.text.split()])
525 stress *= -0.1 * GPa
526 stress = stress.reshape(9)[[0, 4, 8, 5, 2, 1]]
528 dipole = None
529 if lastdipole is not None:
530 dblock = lastdipole.find('v[@name="dipole"]')
531 if dblock is not None:
532 dipole = np.zeros((1, 3), dtype=float)
533 dipole = np.array([float(val) for val in dblock.text.split()])
535 dblock = step.find('dipole/v[@name="dipole"]')
536 if dblock is not None:
537 dipole = np.zeros((1, 3), dtype=float)
538 dipole = np.array([float(val) for val in dblock.text.split()])
540 efermi = step.find('dos/i[@name="efermi"]')
541 if efermi is not None:
542 efermi = float(efermi.text)
544 kpoints = []
545 for ikpt in range(1, len(ibz_kpts) + 1):
546 kblocks = step.findall(
547 'eigenvalues/array/set/set/set[@comment="kpoint %d"]' % ikpt)
548 if kblocks is not None:
549 for spin, kpoint in enumerate(kblocks):
550 eigenvals = kpoint.findall('r')
551 eps_n = np.zeros(len(eigenvals))
552 f_n = np.zeros(len(eigenvals))
553 for j, val in enumerate(eigenvals):
554 val = val.text.split()
555 eps_n[j] = float(val[0])
556 f_n[j] = float(val[1])
557 if len(kblocks) == 1:
558 f_n *= 2
559 kpoints.append(
560 SinglePointKPoint(kpt_weights[ikpt - 1], spin, ikpt,
561 eps_n, f_n))
562 if len(kpoints) == 0:
563 kpoints = None
565 atoms = atoms_init.copy()
566 atoms.set_cell(cell)
567 atoms.set_scaled_positions(scpos)
568 atoms.calc = SinglePointDFTCalculator(atoms,
569 energy=energy,
570 forces=forces,
571 stress=stress,
572 free_energy=free_energy,
573 ibzkpts=ibz_kpts,
574 efermi=efermi,
575 dipole=dipole)
576 atoms.calc.name = 'vasp'
577 atoms.calc.kpts = kpoints
578 atoms.calc.parameters = parameters
579 yield atoms
582@writer
583def write_vasp_xdatcar(fd, images, label=None):
584 """Write VASP MD trajectory (XDATCAR) file
586 Only Vasp 5 format is supported (for consistency with read_vasp_xdatcar)
588 Args:
589 fd (str, fp): Output file
590 images (iterable of Atoms): Atoms images to write. These must have
591 consistent atom order and lattice vectors - this will not be
592 checked.
593 label (str): Text for first line of file. If empty, default to list of
594 elements.
596 """
598 images = iter(images)
599 image = next(images)
601 if not isinstance(image, Atoms):
602 raise TypeError("images should be a sequence of Atoms objects.")
604 symbol_count = _symbol_count_from_symbols(image.get_chemical_symbols())
606 if label is None:
607 label = ' '.join([s for s, _ in symbol_count])
608 fd.write(label + '\n')
610 # Not using lattice constants, set it to 1
611 fd.write(' 1\n')
613 # Lattice vectors; use first image
614 float_string = '{:11.6f}'
615 for row_i in range(3):
616 fd.write(' ')
617 fd.write(' '.join(float_string.format(x) for x in image.cell[row_i]))
618 fd.write('\n')
620 _write_symbol_count(fd, symbol_count)
621 _write_xdatcar_config(fd, image, index=1)
622 for i, image in enumerate(images):
623 # Index is off by 2: 1-indexed file vs 0-indexed Python;
624 # and we already wrote the first block.
625 _write_xdatcar_config(fd, image, i + 2)
628def _write_xdatcar_config(fd, atoms, index):
629 """Write a block of positions for XDATCAR file
631 Args:
632 fd (fd): writeable Python file descriptor
633 atoms (ase.Atoms): Atoms to write
634 index (int): configuration number written to block header
636 """
637 fd.write("Direct configuration={:6d}\n".format(index))
638 float_string = '{:11.8f}'
639 scaled_positions = atoms.get_scaled_positions()
640 for row in scaled_positions:
641 fd.write(' ')
642 fd.write(' '.join([float_string.format(x) for x in row]))
643 fd.write('\n')
646def _symbol_count_from_symbols(symbols):
647 """Reduce list of chemical symbols into compact VASP notation
649 args:
650 symbols (iterable of str)
652 returns:
653 list of pairs [(el1, c1), (el2, c2), ...]
654 """
655 sc = []
656 psym = symbols[0]
657 count = 0
658 for sym in symbols:
659 if sym != psym:
660 sc.append((psym, count))
661 psym = sym
662 count = 1
663 else:
664 count += 1
665 sc.append((psym, count))
666 return sc
669def _write_symbol_count(fd, sc, vasp5=True):
670 """Write the symbols and numbers block for POSCAR or XDATCAR
672 Args:
673 f (fd): Descriptor for writable file
674 sc (list of 2-tuple): list of paired elements and counts
675 vasp5 (bool): if False, omit symbols and only write counts
677 e.g. if sc is [(Sn, 4), (S, 6)] then write::
679 Sn S
680 4 6
682 """
683 if vasp5:
684 for sym, _ in sc:
685 fd.write(' {:3s}'.format(sym))
686 fd.write('\n')
688 for _, count in sc:
689 fd.write(' {:3d}'.format(count))
690 fd.write('\n')
693@writer
694def write_vasp(filename,
695 atoms,
696 label=None,
697 direct=False,
698 sort=None,
699 symbol_count=None,
700 long_format=True,
701 vasp5=True,
702 ignore_constraints=False,
703 wrap=False):
704 """Method to write VASP position (POSCAR/CONTCAR) files.
706 Writes label, scalefactor, unitcell, # of various kinds of atoms,
707 positions in cartesian or scaled coordinates (Direct), and constraints
708 to file. Cartesian coordinates is default and default label is the
709 atomic species, e.g. 'C N H Cu'.
710 """
712 from ase.constraints import FixAtoms, FixScaled, FixedPlane, FixedLine
714 fd = filename # @writer decorator ensures this arg is a file descriptor
716 if isinstance(atoms, (list, tuple)):
717 if len(atoms) > 1:
718 raise RuntimeError('Don\'t know how to save more than ' +
719 'one image to VASP input')
720 else:
721 atoms = atoms[0]
723 # Check lattice vectors are finite
724 if np.any(atoms.cell.cellpar() == 0.):
725 raise RuntimeError(
726 'Lattice vectors must be finite and not coincident. '
727 'At least one lattice length or angle is zero.')
729 # Write atom positions in scaled or cartesian coordinates
730 if direct:
731 coord = atoms.get_scaled_positions(wrap=wrap)
732 else:
733 coord = atoms.get_positions(wrap=wrap)
735 constraints = atoms.constraints and not ignore_constraints
737 if constraints:
738 sflags = np.zeros((len(atoms), 3), dtype=bool)
739 for constr in atoms.constraints:
740 if isinstance(constr, FixScaled):
741 sflags[constr.a] = constr.mask
742 elif isinstance(constr, FixAtoms):
743 sflags[constr.index] = [True, True, True]
744 elif isinstance(constr, FixedPlane):
745 mask = np.all(np.abs(np.cross(constr.dir, atoms.cell)) < 1e-5,
746 axis=1)
747 if sum(mask) != 1:
748 raise RuntimeError(
749 'VASP requires that the direction of FixedPlane '
750 'constraints is parallel with one of the cell axis')
751 sflags[constr.a] = mask
752 elif isinstance(constr, FixedLine):
753 mask = np.all(np.abs(np.cross(constr.dir, atoms.cell)) < 1e-5,
754 axis=1)
755 if sum(mask) != 1:
756 raise RuntimeError(
757 'VASP requires that the direction of FixedLine '
758 'constraints is parallel with one of the cell axis')
759 sflags[constr.a] = ~mask
761 if sort:
762 ind = np.argsort(atoms.get_chemical_symbols())
763 symbols = np.array(atoms.get_chemical_symbols())[ind]
764 coord = coord[ind]
765 if constraints:
766 sflags = sflags[ind]
767 else:
768 symbols = atoms.get_chemical_symbols()
770 # Create a list sc of (symbol, count) pairs
771 if symbol_count:
772 sc = symbol_count
773 else:
774 sc = _symbol_count_from_symbols(symbols)
776 # Create the label
777 if label is None:
778 label = ''
779 for sym, c in sc:
780 label += '%2s ' % sym
781 fd.write(label + '\n')
783 # Write unitcell in real coordinates and adapt to VASP convention
784 # for unit cell
785 # ase Atoms doesn't store the lattice constant separately, so always
786 # write 1.0.
787 fd.write('%19.16f\n' % 1.0)
788 if long_format:
789 latt_form = ' %21.16f'
790 else:
791 latt_form = ' %11.6f'
792 for vec in atoms.get_cell():
793 fd.write(' ')
794 for el in vec:
795 fd.write(latt_form % el)
796 fd.write('\n')
798 # Write out symbols (if VASP 5.x) and counts of atoms
799 _write_symbol_count(fd, sc, vasp5=vasp5)
801 if constraints:
802 fd.write('Selective dynamics\n')
804 if direct:
805 fd.write('Direct\n')
806 else:
807 fd.write('Cartesian\n')
809 if long_format:
810 cform = ' %19.16f'
811 else:
812 cform = ' %9.6f'
813 for iatom, atom in enumerate(coord):
814 for dcoord in atom:
815 fd.write(cform % dcoord)
816 if constraints:
817 for flag in sflags[iatom]:
818 if flag:
819 s = 'F'
820 else:
821 s = 'T'
822 fd.write('%4s' % s)
823 fd.write('\n')