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
|
# fmt: off
import numpy as np
import ase.units
from ase.atoms import Atoms
from ase.calculators.singlepoint import SinglePointCalculator
from ase.data import chemical_symbols
from ase.io import read
def parse_geometry(filename):
'''Read atoms geometry from ACE-Molecule log file and put it to self.data.
Parameters
==========
filename: ACE-Molecule log file.
Returns
=======
Dictionary of parsed geometry data.
retval["Atomic_numbers"]: list of atomic numbers
retval["Positions"]: list of [x, y, z] coordinates for each atoms.
'''
with open(filename) as fd:
lines = fd.readlines()
start_line = 0
end_line = 0
for i, line in enumerate(lines):
if line == '==================== Atoms =====================\n':
start_line = i
if start_line != 0 and len(line.split('=')) > 3:
end_line = i
if start_line != end_line:
break
atoms = []
positions = []
for i in range(start_line + 1, end_line):
atomic_number = lines[i].split()[0]
atoms.append(str(chemical_symbols[int(atomic_number)]))
xyz = [float(n) for n in lines[i].split()[1:4]]
positions.append(xyz)
return {"Atomic_numbers": atoms, "Positions": positions}
def read_acemolecule_out(filename):
'''Interface to ACEMoleculeReader, return values for corresponding quantity
Parameters
==========
filename: ACE-Molecule log file.
quantity: One of atoms, energy, forces, excitation-energy.
Returns
=======
- quantity = 'excitation-energy':
returns None. This is placeholder function to run TDDFT calculations
without IndexError.
- quantity = 'energy':
returns energy as float value.
- quantity = 'forces':
returns force of each atoms as numpy array of shape (natoms, 3).
- quantity = 'atoms':
returns ASE atoms object.
'''
data = parse_geometry(filename)
atom_symbol = np.array(data["Atomic_numbers"])
positions = np.array(data["Positions"])
atoms = Atoms(atom_symbol, positions=positions)
energy = None
forces = None
excitation_energy = None
# results = {}
# if len(results)<1:
with open(filename) as fd:
lines = fd.readlines()
for i in range(len(lines) - 1, 1, -1):
line = lines[i].split()
if len(line) > 2:
if line[0] == 'Total' and line[1] == 'energy':
energy = float(line[3])
break
# energy must be modified, hartree to eV
energy *= ase.units.Hartree
forces = []
for i in range(len(lines) - 1, 1, -1):
if '!============================' in lines[i]:
endline_num = i
if '! Force:: List of total force in atomic unit' in lines[i]:
startline_num = i + 2
for j in range(startline_num, endline_num):
forces.append(lines[j].split()[3:6])
convert = ase.units.Hartree / ase.units.Bohr
forces = np.array(forces, dtype=float) * convert
break
if len(forces) <= 0:
forces = None
# Set calculator to
calc = SinglePointCalculator(atoms, energy=energy, forces=forces)
atoms.calc = calc
results = {}
results['energy'] = energy
results['atoms'] = atoms
results['forces'] = forces
results['excitation-energy'] = excitation_energy
return results
def read_acemolecule_input(filename):
'''Reads a ACE-Molecule input file
Parameters
==========
filename: ACE-Molecule input file name
Returns
=======
ASE atoms object containing geometry only.
'''
with open(filename) as fd:
for line in fd:
if len(line.split('GeometryFilename')) > 1:
geometryfile = line.split()[1]
break
atoms = read(geometryfile, format='xyz')
return atoms
|
