Source code for

from typing import Sequence, Union

import numpy as np

from ase.atoms import Atoms
from ase.utils import reader, writer

[docs]@reader def read_dftb(fd): """Method to read coordinates from the Geometry section of a DFTB+ input file (typically called "dftb_in.hsd"). As described in the DFTB+ manual, this section can be in a number of different formats. This reader supports the GEN format and the so-called "explicit" format. The "explicit" format is unique to DFTB+ input files. The GEN format can also be used in a stand-alone fashion, as coordinate files with a `.gen` extension. Reading and writing such files is implemented in ``. """ lines = fd.readlines() atoms_pos = [] atom_symbols = [] type_names = [] my_pbc = False fractional = False mycell = [] for iline, line in enumerate(lines): if line.strip().startswith('#'): pass elif 'genformat' in line.lower(): natoms = int(lines[iline + 1].split()[0]) if lines[iline + 1].split()[1].lower() == 's': my_pbc = True elif lines[iline + 1].split()[1].lower() == 'f': my_pbc = True fractional = True symbols = lines[iline + 2].split() for i in range(natoms): index = iline + 3 + i aindex = int(lines[index].split()[1]) - 1 atom_symbols.append(symbols[aindex]) position = [float(p) for p in lines[index].split()[2:]] atoms_pos.append(position) if my_pbc: for i in range(3): index = iline + 4 + natoms + i cell = [float(c) for c in lines[index].split()] mycell.append(cell) else: if 'TypeNames' in line: col = line.split() for i in range(3, len(col) - 1): type_names.append(col[i].strip("\"")) elif 'Periodic' in line: if 'Yes' in line: my_pbc = True elif 'LatticeVectors' in line: for imycell in range(3): extraline = lines[iline + imycell + 1] cols = extraline.split() mycell.append( [float(cols[0]), float(cols[1]), float(cols[2])]) else: pass if not my_pbc: mycell = [0.] * 3 start_reading_coords = False stop_reading_coords = False for line in lines: if line.strip().startswith('#'): pass else: if 'TypesAndCoordinates' in line: start_reading_coords = True if start_reading_coords: if '}' in line: stop_reading_coords = True if (start_reading_coords and not stop_reading_coords and 'TypesAndCoordinates' not in line): typeindexstr, xxx, yyy, zzz = line.split()[:4] typeindex = int(typeindexstr) symbol = type_names[typeindex - 1] atom_symbols.append(symbol) atoms_pos.append([float(xxx), float(yyy), float(zzz)]) if fractional: atoms = Atoms(scaled_positions=atoms_pos, symbols=atom_symbols, cell=mycell, pbc=my_pbc) elif not fractional: atoms = Atoms(positions=atoms_pos, symbols=atom_symbols, cell=mycell, pbc=my_pbc) return atoms
def read_dftb_velocities(atoms, filename): """Method to read velocities (AA/ps) from DFTB+ output file """ from ase.units import second # AA/ps -> ase units AngdivPs2ASE = 1.0 / (1e-12 * second) with open(filename) as fd: lines = fd.readlines() # remove empty lines lines_ok = [] for line in lines: if line.rstrip(): lines_ok.append(line) velocities = [] natoms = len(atoms) last_lines = lines_ok[-natoms:] for iline, line in enumerate(last_lines): inp = line.split() velocities.append([float(inp[5]) * AngdivPs2ASE, float(inp[6]) * AngdivPs2ASE, float(inp[7]) * AngdivPs2ASE]) atoms.set_velocities(velocities) return atoms @reader def read_dftb_lattice(fileobj, images=None): """Read lattice vectors from MD and return them as a list. If a molecules are parsed add them there.""" if images is not None: append = True if hasattr(images, 'get_positions'): images = [images] else: append = False lattices = [] for line in fileobj: if 'Lattice vectors' in line: vec = [] for _ in range(3): line = fileobj.readline().split() try: line = [float(x) for x in line] except ValueError: raise ValueError('Lattice vector elements should be of ' 'type float.') vec.extend(line) lattices.append(np.array(vec).reshape((3, 3))) if append: if len(images) != len(lattices): raise ValueError('Length of images given does not match number of ' 'cell vectors found') for i, atoms in enumerate(images): atoms.set_cell(lattices[i]) # DFTB+ only supports 3D PBC atoms.set_pbc(True) return None else: return lattices
[docs]@writer def write_dftb( fileobj, images: Union[Atoms, Sequence[Atoms]], fractional: bool = False, ): """Write structure in GEN format (refer to DFTB+ manual). Multiple snapshots are not allowed. """ from import write_gen write_gen(fileobj, images, fractional=fractional)
def write_dftb_velocities(atoms, filename): """Method to write velocities (in atomic units) from ASE to a file to be read by dftb+ """ from ase.units import AUT, Bohr # ase units -> atomic units ASE2au = Bohr / AUT with open(filename, 'w') as fd: velocities = atoms.get_velocities() for velocity in velocities: fd.write(' %19.16f %19.16f %19.16f \n' % (velocity[0] / ASE2au, velocity[1] / ASE2au, velocity[2] / ASE2au))