|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +"""Command line utilities to create and inspect `StructureData` nodes from CP2K input files.""" |
| 3 | + |
| 4 | +from fractions import Fraction |
| 5 | + |
| 6 | +import click |
| 7 | +import numpy as np |
| 8 | + |
| 9 | +from aiida.cmdline.params import options |
| 10 | +from aiida.cmdline.utils import decorators, echo |
| 11 | + |
| 12 | +from . import cmd_data |
| 13 | + |
| 14 | + |
| 15 | +@cmd_data.group('structure') |
| 16 | +def cmd_structure(): |
| 17 | + """Commands to create and inspect `StructureData` nodes from CP2K input.""" |
| 18 | + |
| 19 | + |
| 20 | +@cmd_structure.command('import') |
| 21 | +@click.argument('filename', type=click.File('r')) |
| 22 | +@options.DRY_RUN() |
| 23 | +@decorators.with_dbenv() |
| 24 | +def cmd_import(filename, dry_run): |
| 25 | + """Import a `StructureData` from a CP2K input file.""" |
| 26 | + |
| 27 | + # pylint: disable=import-outside-toplevel,invalid-name,too-many-locals,too-many-statements,too-many-branches |
| 28 | + |
| 29 | + from cp2k_input_tools.parser import CP2KInputParser |
| 30 | + from aiida.orm.nodes.data.structure import StructureData, Kind, Site, symop_fract_from_ortho |
| 31 | + |
| 32 | + parser = CP2KInputParser() |
| 33 | + tree = parser.parse(filename) |
| 34 | + |
| 35 | + for force_eval in tree['+force_eval']: |
| 36 | + try: |
| 37 | + cell = force_eval['+subsys']['+cell'] |
| 38 | + kinds = force_eval['+subsys']['+kind'] |
| 39 | + coord = force_eval['+subsys']['+coord'] |
| 40 | + break # for now grab the first &COORD found |
| 41 | + except KeyError: |
| 42 | + continue |
| 43 | + else: |
| 44 | + echo.echo_error('no STRUCTURE, CELL, or KIND found in the given input file') |
| 45 | + |
| 46 | + structure = StructureData() |
| 47 | + |
| 48 | + # CP2K can get its cell information in two ways: |
| 49 | + # - A, B, C: cell vectors |
| 50 | + # - ABC: scaling of cell vectors, ALPHA_BETA_GAMMA: angles between the cell vectors |
| 51 | + # We'll parse either form |
| 52 | + if 'a' in cell: |
| 53 | + unit_cell = np.array([cell['a'], cell['b'], cell['c']]) # unit vectors given |
| 54 | + |
| 55 | + norms = [np.linalg.norm(v) for v in unit_cell] |
| 56 | + angles = [] |
| 57 | + |
| 58 | + for i in range(3): |
| 59 | + j = i - 1 |
| 60 | + k = i - 2 |
| 61 | + |
| 62 | + ll = norms[j] * norms[k] |
| 63 | + if ll > 1e-16: |
| 64 | + angle = np.arccos(np.dot(unit_cell[j], unit_cell[k]) / ll) |
| 65 | + else: |
| 66 | + angle = np.pi / 2 |
| 67 | + |
| 68 | + angles.append(angle) |
| 69 | + |
| 70 | + tmat = symop_fract_from_ortho(norms + angles) |
| 71 | + |
| 72 | + elif 'abc' in cell: |
| 73 | + # length of unit vectors given |
| 74 | + if 'alpha_beta_gamma' in cell: |
| 75 | + # if we also have the angles, construct the cell |
| 76 | + alpha, beta, gamma = cell.pop('alpha_beta_gamma') |
| 77 | + |
| 78 | + cos_alpha = np.cos(alpha) |
| 79 | + cos_beta = np.cos(beta) |
| 80 | + cos_gamma = np.cos(gamma) |
| 81 | + sin_gamma = np.sin(gamma) |
| 82 | + |
| 83 | + unit_cell = np.array([ |
| 84 | + [1., 0., 0.], |
| 85 | + [cos_gamma, sin_gamma, 0.], |
| 86 | + [ |
| 87 | + cos_beta, |
| 88 | + (cos_alpha - cos_gamma * cos_beta) / sin_gamma, |
| 89 | + np.sqrt(1. - cos_beta**2 - ((cos_alpha - cos_gamma * cos_beta) / sin_gamma)**2), |
| 90 | + ], |
| 91 | + ]) |
| 92 | + else: |
| 93 | + unit_cell = np.eye(3) |
| 94 | + alpha = beta = gamma = np.pi / 2 |
| 95 | + |
| 96 | + a, b, c = cell['abc'] # unpack and scale |
| 97 | + unit_cell[0, :] *= a |
| 98 | + unit_cell[1, :] *= b |
| 99 | + unit_cell[2, :] *= c |
| 100 | + |
| 101 | + tmat = symop_fract_from_ortho([a, b, c, alpha, beta, gamma]) |
| 102 | + else: |
| 103 | + echo.echo_error('incomplete &CELL section') |
| 104 | + |
| 105 | + if not coord.get('scaled', False): |
| 106 | + tmat = np.eye(3) |
| 107 | + |
| 108 | + structure.set_attribute('cell', unit_cell) |
| 109 | + |
| 110 | + for kind in kinds: |
| 111 | + name = kind['_'] |
| 112 | + element = kind.get('element', name) # TODO: support getting the element from e.g. 'Fe1' |
| 113 | + structure.append_kind(Kind(name=name, symbols=element)) |
| 114 | + |
| 115 | + if coord.get('units', 'angstrom').lower() != 'angstrom': |
| 116 | + echo.echo_error('unit conversion for coordinations is not (yet) supported') |
| 117 | + |
| 118 | + for coordline in coord['*']: |
| 119 | + # coordinates are a series of strings according to the CP2K schema |
| 120 | + fields = coordline.split() |
| 121 | + |
| 122 | + name = fields[0] |
| 123 | + # positions can be fractions AND they may be scaled |
| 124 | + position = tmat @ np.array([float(Fraction(f)) for f in fields[1:4]]) |
| 125 | + |
| 126 | + structure.append_site(Site(kind_name=name, position=position)) |
| 127 | + |
| 128 | + formula = structure.get_formula() |
| 129 | + |
| 130 | + if dry_run: |
| 131 | + echo.echo_success('parsed structure with formula {}'.format(formula)) |
| 132 | + else: |
| 133 | + structure.store() |
| 134 | + echo.echo_success('parsed and stored StructureData<{}> with formula {}'.format(structure.pk, formula)) |
0 commit comments