File: gen.py

package info (click to toggle)
python-ase 3.24.0-1
  • links: PTS, VCS
  • area: main
  • in suites: trixie
  • size: 15,448 kB
  • sloc: python: 144,945; xml: 2,728; makefile: 113; javascript: 47
file content (137 lines) | stat: -rw-r--r-- 4,608 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
"""Extension to ASE: read and write structures in GEN format

Refer to DFTB+ manual for GEN format description.

Note: GEN format only supports single snapshot.
"""
from typing import Dict, Sequence, Union

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


@reader
def read_gen(fileobj):
    """Read structure in GEN format (refer to DFTB+ manual).
       Multiple snapshot are not allowed. """
    image = Atoms()
    lines = fileobj.readlines()
    line = lines[0].split()
    natoms = int(line[0])
    pb_flag = line[1]
    if line[1] not in ['C', 'F', 'S']:
        if line[1] == 'H':
            raise OSError('Error in line #1: H (Helical) is valid but not '
                          'supported. Only C (Cluster), S (Supercell) '
                          'or F (Fraction) are supported options')
        else:
            raise OSError('Error in line #1: only C (Cluster), S (Supercell) '
                          'or F (Fraction) are supported options')

    # Read atomic symbols
    line = lines[1].split()
    symboldict = {symbolid: symb for symbolid, symb in enumerate(line, start=1)}
    # Read atoms (GEN format supports only single snapshot)
    del lines[:2]
    positions = []
    symbols = []
    for line in lines[:natoms]:
        _dummy, symbolid, x, y, z = line.split()[:5]
        symbols.append(symboldict[int(symbolid)])
        positions.append([float(x), float(y), float(z)])
    image = Atoms(symbols=symbols, positions=positions)
    del lines[:natoms]

    # If Supercell, parse periodic vectors.
    # If Fraction, translate into Supercell.
    if pb_flag == 'C':
        return image
    else:
        # Dummy line: line after atom positions is not uniquely defined
        # in gen implementations, and not necessary in DFTB package
        del lines[:1]
        image.set_pbc([True, True, True])
        p = []
        for i in range(3):
            x, y, z = lines[i].split()[:3]
            p.append([float(x), float(y), float(z)])
        image.set_cell([(p[0][0], p[0][1], p[0][2]),
                        (p[1][0], p[1][1], p[1][2]),
                        (p[2][0], p[2][1], p[2][2])])
        if pb_flag == 'F':
            frac_positions = image.get_positions()
            image.set_scaled_positions(frac_positions)
        return image


@writer
def write_gen(
    fileobj,
    images: Union[Atoms, Sequence[Atoms]],
    fractional: bool = False,
):
    """Write structure in GEN format (refer to DFTB+ manual).
       Multiple snapshots are not allowed. """
    if isinstance(images, (list, tuple)):
        # GEN format doesn't support multiple snapshots
        if len(images) != 1:
            raise ValueError(
                '"images" contains more than one structure. '
                'GEN format supports only single snapshot output.'
            )
        atoms = images[0]
    else:
        atoms = images

    symbols = atoms.get_chemical_symbols()

    # Define a dictionary with symbols-id
    symboldict: Dict[str, int] = {}
    for sym in symbols:
        if sym not in symboldict:
            symboldict[sym] = len(symboldict) + 1
    # An ordered symbol list is needed as ordered dictionary
    # is just available in python 2.7
    orderedsymbols = list(['null'] * len(symboldict.keys()))
    for sym, num in symboldict.items():
        orderedsymbols[num - 1] = sym

    # Check whether the structure is periodic
    # GEN cannot describe periodicity in one or two direction,
    # a periodic structure is considered periodic in all the
    # directions. If your structure is not periodical in all
    # the directions, be sure you have set big periodicity
    # vectors in the non-periodic directions
    if fractional:
        pb_flag = 'F'
    elif atoms.pbc.any():
        pb_flag = 'S'
    else:
        pb_flag = 'C'

    natoms = len(symbols)
    ind = 0

    fileobj.write(f'{natoms:d}  {pb_flag:<5s}\n')
    for sym in orderedsymbols:
        fileobj.write(f'{sym:<5s}')
    fileobj.write('\n')

    if fractional:
        coords = atoms.get_scaled_positions(wrap=False)
    else:
        coords = atoms.get_positions(wrap=False)

    for sym, (x, y, z) in zip(symbols, coords):
        ind += 1
        symbolid = symboldict[sym]
        fileobj.write(
            f'{ind:-6d} {symbolid:d} {x:22.15f} {y:22.15f} {z:22.15f}\n')

    if atoms.pbc.any() or fractional:
        fileobj.write(f'{0.0:22.15f} {0.0:22.15f} {0.0:22.15f} \n')
        cell = atoms.get_cell()
        for i in range(3):
            for j in range(3):
                fileobj.write(f'{cell[i, j]:22.15f} ')
            fileobj.write('\n')