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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
|
"""This module defines an ASE interface to MOPAC.
Set $ASE_MOPAC_COMMAND to something like::
LD_LIBRARY_PATH=/path/to/lib/ \
MOPAC_LICENSE=/path/to/license \
/path/to/MOPAC2012.exe PREFIX.mop 2> /dev/null
"""
import os
import numpy as np
from ase import Atoms
from ase.calculators.calculator import FileIOCalculator, ReadError, Parameters
from ase.units import kcal, mol, Debye
class MOPAC(FileIOCalculator):
implemented_properties = ['energy', 'forces', 'dipole', 'magmom']
command = 'mopac PREFIX.mop 2> /dev/null'
discard_results_on_any_change = True
default_parameters = dict(
method='PM7',
task='1SCF GRADIENTS',
charge=0,
relscf=0.0001)
methods = ['AM1', 'MNDO', 'MNDOD', 'PM3', 'PM6', 'PM6-D3', 'PM6-DH+',
'PM6-DH2', 'PM6-DH2X', 'PM6-D3H4', 'PM6-D3H4X', 'PMEP', 'PM7',
'PM7-TS', 'RM1']
def __init__(self, restart=None,
ignore_bad_restart_file=FileIOCalculator._deprecated,
label='mopac', atoms=None, **kwargs):
"""Construct MOPAC-calculator object.
Parameters:
label: str
Prefix for filenames (label.mop, label.out, ...)
Examples:
Use default values to do a single SCF calculation and print
the forces (task='1SCF GRADIENTS'):
>>> from ase.build import molecule
>>> from ase.calculators.mopac import MOPAC
>>> atoms = molecule('O2')
>>> atoms.calc = MOPAC(label='O2')
>>> atoms.get_potential_energy()
>>> eigs = atoms.calc.get_eigenvalues()
>>> somos = atoms.calc.get_somo_levels()
>>> homo, lumo = atoms.calc.get_homo_lumo_levels()
Use the internal geometry optimization of Mopac:
>>> atoms = molecule('H2')
>>> atoms.calc = MOPAC(label='H2', task='GRADIENTS')
>>> atoms.get_potential_energy()
Read in and start from output file:
>>> atoms = MOPAC.read_atoms('H2')
>>> atoms.calc.get_homo_lumo_levels()
"""
FileIOCalculator.__init__(self, restart, ignore_bad_restart_file,
label, atoms, **kwargs)
def write_input(self, atoms, properties=None, system_changes=None):
FileIOCalculator.write_input(self, atoms, properties, system_changes)
p = self.parameters
# Build string to hold .mop input file:
s = p.method + ' ' + p.task + ' '
if p.relscf:
s += 'RELSCF={0} '.format(p.relscf)
# Write charge:
if p.charge:
charge = p.charge
else:
charge = atoms.get_initial_charges().sum()
if charge != 0:
s += 'CHARGE={0} '.format(int(round(charge)))
magmom = int(round(abs(atoms.get_initial_magnetic_moments().sum())))
if magmom:
s += (['DOUBLET', 'TRIPLET', 'QUARTET', 'QUINTET'][magmom - 1] +
' UHF ')
s += '\nTitle: ASE calculation\n\n'
# Write coordinates:
for xyz, symbol in zip(atoms.positions, atoms.get_chemical_symbols()):
s += ' {0:2} {1} 1 {2} 1 {3} 1\n'.format(symbol, *xyz)
for v, p in zip(atoms.cell, atoms.pbc):
if p:
s += 'Tv {0} {1} {2}\n'.format(*v)
with open(self.label + '.mop', 'w') as f:
f.write(s)
def get_spin_polarized(self):
return self.nspins == 2
def get_index(self, lines, pattern):
for i, line in enumerate(lines):
if line.find(pattern) != -1:
return i
def read(self, label):
FileIOCalculator.read(self, label)
if not os.path.isfile(self.label + '.out'):
raise ReadError
with open(self.label + '.out') as f:
lines = f.readlines()
self.parameters = Parameters(task='', method='')
p = self.parameters
parm_line = self.read_parameters_from_file(lines)
for keyword in parm_line.split():
if 'RELSCF' in keyword:
p.relscf = float(keyword.split('=')[-1])
elif keyword in self.methods:
p.method = keyword
else:
p.task += keyword + ' '
p.task.rstrip()
self.atoms = self.read_atoms_from_file(lines)
self.read_results()
def read_atoms_from_file(self, lines):
"""Read the Atoms from the output file stored as list of str in lines.
Parameters:
lines: list of str
"""
# first try to read from final point (last image)
i = self.get_index(lines, 'FINAL POINT AND DERIVATIVES')
if i is None: # XXX should we read it from the input file?
assert 0, 'Not implemented'
lines1 = lines[i:]
i = self.get_index(lines1, 'CARTESIAN COORDINATES')
j = i + 2
symbols = []
positions = []
while not lines1[j].isspace(): # continue until we hit a blank line
l = lines1[j].split()
symbols.append(l[1])
positions.append([float(c) for c in l[2: 2 + 3]])
j += 1
return Atoms(symbols=symbols, positions=positions)
def read_parameters_from_file(self, lines):
"""Find and return the line that defines a Mopac calculation
Parameters:
lines: list of str
"""
for i, line in enumerate(lines):
if line.find('CALCULATION DONE:') != -1:
break
lines1 = lines[i:]
for i, line in enumerate(lines1):
if line.find('****') != -1:
return lines1[i + 1]
def read_results(self):
"""Read the results, such as energy, forces, eigenvalues, etc.
"""
FileIOCalculator.read(self, self.label)
if not os.path.isfile(self.label + '.out'):
raise ReadError
with open(self.label + '.out') as f:
lines = f.readlines()
for i, line in enumerate(lines):
if line.find('TOTAL ENERGY') != -1:
self.results['energy'] = float(line.split()[3])
elif line.find('FINAL HEAT OF FORMATION') != -1:
self.final_hof = float(line.split()[5]) * kcal / mol
elif line.find('NO. OF FILLED LEVELS') != -1:
self.nspins = 1
self.no_occ_levels = int(line.split()[-1])
elif line.find('NO. OF ALPHA ELECTRON') != -1:
self.nspins = 2
self.no_alpha_electrons = int(line.split()[-1])
self.no_beta_electrons = int(lines[i+1].split()[-1])
self.results['magmom'] = abs(self.no_alpha_electrons -
self.no_beta_electrons)
elif line.find('FINAL POINT AND DERIVATIVES') != -1:
forces = [-float(line.split()[6])
for line in lines[i + 3:i + 3 + 3 * len(self.atoms)]]
self.results['forces'] = np.array(
forces).reshape((-1, 3)) * kcal / mol
elif line.find('EIGENVALUES') != -1:
if line.find('ALPHA') != -1:
j = i + 1
eigs_alpha = []
while not lines[j].isspace():
eigs_alpha += [float(eps) for eps in lines[j].split()]
j += 1
elif line.find('BETA') != -1:
j = i + 1
eigs_beta = []
while not lines[j].isspace():
eigs_beta += [float(eps) for eps in lines[j].split()]
j += 1
eigs = np.array([eigs_alpha, eigs_beta]).reshape(2, 1, -1)
self.eigenvalues = eigs
else:
eigs = []
j = i + 1
while not lines[j].isspace():
eigs += [float(e) for e in lines[j].split()]
j += 1
self.eigenvalues = np.array(eigs).reshape(1, 1, -1)
elif line.find('DIPOLE ') != -1:
self.results['dipole'] = np.array(
lines[i + 3].split()[1:1 + 3], float) * Debye
def get_eigenvalues(self, kpt=0, spin=0):
return self.eigenvalues[spin, kpt]
def get_homo_lumo_levels(self):
eigs = self.eigenvalues
if self.nspins == 1:
nocc = self.no_occ_levels
return np.array([eigs[0, 0, nocc - 1], eigs[0, 0, nocc]])
else:
na = self.no_alpha_electrons
nb = self.no_beta_electrons
if na == 0:
return None, self.eigenvalues[1, 0, nb - 1]
elif nb == 0:
return self.eigenvalues[0, 0, na - 1], None
else:
eah, eal = eigs[0, 0, na - 1: na + 1]
ebh, ebl = eigs[1, 0, nb - 1: nb + 1]
return np.array([max(eah, ebh), min(eal, ebl)])
def get_somo_levels(self):
assert self.nspins == 2
na, nb = self.no_alpha_electrons, self.no_beta_electrons
if na == 0:
return None, self.eigenvalues[1, 0, nb - 1]
elif nb == 0:
return self.eigenvalues[0, 0, na - 1], None
else:
return np.array([self.eigenvalues[0, 0, na - 1],
self.eigenvalues[1, 0, nb - 1]])
def get_final_heat_of_formation(self):
"""Final heat of formation as reported in the Mopac output file
"""
return self.final_hof
|
