Source code for ase.io.nwchem.nwreader_in

import re

import numpy as np

from ase import Atoms
from ase.geometry import cellpar_to_cell

from .parser import _define_pattern

# Geometry block parser
_geom = _define_pattern(
    r'^[ \t]*geometry[ \t\S]*\n'
    r'((?:^[ \t]*[\S]+[ \t\S]*\n)+)'
    r'^[ \t]*end\n\n',
    """\
geometry units angstrom nocenter noautosym noautoz
  system crystal units angstrom
    lattice_vectors
      4.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
      0.0000000000000000e+00 5.5264780000000000e+00 0.0000000000000000e+00
      0.0000000000000000e+00 0.0000000000000000e+00 4.5963089999999998e+00
  end
   O 5.0000000000000000e-01 5.0000000000000011e-01 5.6486824536818558e-01
   H 5.0000000000000000e-01 6.3810586054988372e-01 4.3513175463181430e-01
   H 5.0000000000000000e-01 3.6189413945011639e-01 4.3513175463181430e-01
end

""", re.M)

# Finds crystal specification
_crystal = _define_pattern(
    r'^[ \t]*system crystal[ \t\S]*\n'
    r'((?:[ \t]*[\S]+[ \t\S]*\n)+?)'
    r'^[ \t]*end[ \t]*\n',
    """\
  system crystal units angstrom
    lattice_vectors
      4.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
      0.0000000000000000e+00 5.5264780000000000e+00 0.0000000000000000e+00
      0.0000000000000000e+00 0.0000000000000000e+00 4.5963089999999998e+00
  end
""", re.M)

# Finds 3d-periodic unit cell
_cell_3d = _define_pattern(
    r'^[ \t]*lattice_vectors[ \t]*\n'
    r'^((?:(?:[ \t]+[\S]+){3}\n){3})',
    """\
    lattice_vectors
      4.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
      0.0000000000000000e+00 5.5264780000000000e+00 0.0000000000000000e+00
      0.0000000000000000e+00 0.0000000000000000e+00 4.5963089999999998e+00
""", re.M)

# Extracts chemical species from a geometry block
_species = _define_pattern(
    r'^[ \t]*[A-Z][a-z]?(?:[ \t]+[\S]+){3}\n',
    "   O 0.0 0.0 0.0\n", re.M,
)


[docs]def read_nwchem_in(fobj, index=-1): text = ''.join(fobj.readlines()) atomslist = [] for match in _geom.findall(text): symbols = [] positions = [] for atom in _species.findall(match): atom = atom.split() symbols.append(atom[0]) positions.append([float(x) for x in atom[1:]]) positions = np.array(positions) atoms = Atoms(symbols) cell, pbc = _get_cell(text) pos = np.zeros_like(positions) for dim, ipbc in enumerate(pbc): if ipbc: pos += np.outer(positions[:, dim], cell[dim, :]) else: pos[:, dim] = positions[:, dim] atoms.set_cell(cell) atoms.pbc = pbc atoms.set_positions(pos) atomslist.append(atoms) return atomslist[index]
def _get_cell(text): # first check whether there is a lattice definition cell = np.zeros((3, 3)) lattice = _cell_3d.findall(text) if lattice: pbc = [True, True, True] for i, row in enumerate(lattice[0].strip().split('\n')): cell[i] = [float(x) for x in row.split()] return cell, pbc pbc = [False, False, False] lengths = [None, None, None] angles = [None, None, None] for row in text.strip().split('\n'): row = row.strip().lower() for dim, vecname in enumerate(['a', 'b', 'c']): if row.startswith(f'lat_{vecname}'): pbc[dim] = True lengths[dim] = float(row.split()[1]) for i, angle in enumerate(['alpha', 'beta', 'gamma']): if row.startswith(angle): angles[i] = float(row.split()[1]) if not np.any(pbc): return None, pbc for i in range(3): a, b, c = np.roll(np.array([0, 1, 2]), i) if pbc[a] and pbc[b]: assert angles[c] is not None if angles[c] is not None: assert pbc[a] and pbc[b] # The easiest case: all three lattice vectors and angles are specified if np.all(pbc): return cellpar_to_cell(lengths + angles), pbc # Next easiest case: exactly one lattice vector has been specified if np.sum(pbc) == 1: dim = np.argmax(pbc) cell[dim, dim] = lengths[dim] return cell, pbc # Hardest case: two lattice vectors are specified. dim1, dim2 = (dim for dim, ipbc in enumerate(pbc) if ipbc) angledim = np.argmin(pbc) cell[dim1, dim1] = lengths[dim1] cell[dim2, dim2] = lengths[dim2] * np.sin(angles[angledim]) cell[dim2, dim1] = lengths[dim2] * np.cos(angles[angledim]) return cell, pbc