Source code for ase.io.dftb

import numpy as np
from ase.atoms import Atoms


[docs]def read_dftb(filename='dftb_in.hsd'): """Method to read coordinates form DFTB+ input file dftb_in.hsd additionally read information about fixed atoms and periodic boundary condition """ with open(filename, 'r') as myfile: lines = myfile.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 = [1.0, 1.0, 1.0] 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='geo_end.xyz'): """Method to read velocities (AA/ps) from DFTB+ output file geo_end.xyz """ from ase.units import second # AA/ps -> ase units AngdivPs2ASE = 1.0 / (1e-12 * second) myfile = open(filename) lines = myfile.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 def read_dftb_lattice(fileobj='md.out', images=None): """Read lattice vectors from MD and return them as a list. If a molecules are parsed add them there.""" if isinstance(fileobj, str): fileobj = open(fileobj) if images is not None: append = True if hasattr(images, 'get_positions'): images = [images] else: append = False fileobj.seek(0) lattices = [] for line in fileobj: if 'Lattice vectors' in line: vec = [] for i in range(3): # DFTB+ only supports 3D PBC 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 else: return lattices def write_dftb_velocities(atoms, filename='velocities.txt'): """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 if isinstance(filename, str): myfile = open(filename, 'w') else: # Assume it's a 'file-like object' myfile = filename velocities = atoms.get_velocities() for velocity in velocities: myfile.write(' %19.16f %19.16f %19.16f \n' % (velocity[0] / ASE2au, velocity[1] / ASE2au, velocity[2] / ASE2au)) return
[docs]def write_dftb(filename, atoms): """Method to write atom structure in DFTB+ format (gen format) """ # sort atoms.set_masses() masses = atoms.get_masses() indexes = np.argsort(masses) atomsnew = Atoms() for i in indexes: atomsnew = atomsnew + atoms[i] if isinstance(filename, str): myfile = open(filename, 'w') else: # Assume it's a 'file-like object' myfile = filename ispbc = atoms.get_pbc() box = atoms.get_cell() if (any(ispbc)): myfile.write('%8d %2s \n' % (len(atoms), 'S')) else: myfile.write('%8d %2s \n' % (len(atoms), 'C')) chemsym = atomsnew.get_chemical_symbols() allchem = '' for i in chemsym: if i not in allchem: allchem = allchem + i + ' ' myfile.write(allchem + ' \n') coords = atomsnew.get_positions() itype = 1 for iatom, coord in enumerate(coords): if iatom > 0: if chemsym[iatom] != chemsym[iatom - 1]: itype = itype + 1 myfile.write('%5i%5i %19.16f %19.16f %19.16f \n' % (iatom + 1, itype, coords[iatom][0], coords[iatom][1], coords[iatom][2])) # write box if (any(ispbc)): # dftb dummy myfile.write(' %19.16f %19.16f %19.16f \n' % (0, 0, 0)) myfile.write(' %19.16f %19.16f %19.16f \n' % (box[0][0], box[0][1], box[0][2])) myfile.write(' %19.16f %19.16f %19.16f \n' % (box[1][0], box[1][1], box[1][2])) myfile.write(' %19.16f %19.16f %19.16f \n' % (box[2][0], box[2][1], box[2][2])) if isinstance(filename, str): myfile.close()