Coverage for ase / calculators / mopac.py: 85.56%

1# fmt: off 

2 

3"""This module defines an ASE interface to MOPAC. 

4 

5Set $ASE_MOPAC_COMMAND to something like:: 

6 

7 LD_LIBRARY_PATH=/path/to/lib/ \ 

8 MOPAC_LICENSE=/path/to/license \ 

9 /path/to/MOPAC2012.exe PREFIX.mop 2> /dev/null 

10 

11""" 

12import os 

13import re 

14from collections.abc import Sequence 

15from warnings import warn 

16 

17import numpy as np 

18from packaging import version 

19 

20from ase import Atoms 

21from ase.calculators.calculator import FileIOCalculator, Parameters, ReadError 

22from ase.units import Debye, kcal, mol 

23 

24 

25def get_version_number(lines: Sequence[str]): 

26 pattern1 = r'MOPAC\s+v(\S+)' 

27 pattern2 = r'MOPAC2016, Version:\s+([^,]+)' 

28 

29 for line in lines[:10]: 

30 match = re.search(pattern1, line) or re.search(pattern2, line) 

31 if match: 

32 return match.group(1) 

33 raise ValueError('Version number was not found in MOPAC output') 

34 

35 

36class MOPAC(FileIOCalculator): 

37 implemented_properties = ['energy', 'forces', 'dipole', 

38 'magmom', 'free_energy'] 

39 _legacy_default_command = 'mopac PREFIX.mop 2> /dev/null' 

40 discard_results_on_any_change = True 

41 

42 default_parameters = dict( 

43 method='PM7', 

44 task='1SCF GRADIENTS', 

45 charge=None, 

46 relscf=0.0001) 

47 

48 methods = ['AM1', 'MNDO', 'MNDOD', 'PM3', 'PM6', 'PM6-D3', 'PM6-DH+', 

49 'PM6-DH2', 'PM6-DH2X', 'PM6-D3H4', 'PM6-D3H4X', 'PMEP', 'PM7', 

50 'PM7-TS', 'RM1'] 

51 

52 fileio_rules = FileIOCalculator.ruleset( 

53 extend_argv=['{prefix}.mop'], 

54 stdout_name='{prefix}.out') 

55 

56 def __init__(self, restart=None, 

57 ignore_bad_restart_file=FileIOCalculator._deprecated, 

58 label='mopac', atoms=None, **kwargs): 

59 """Construct MOPAC-calculator object. 

60 

61 Parameters 

62 ---------- 

63 

64 label: str 

65 Prefix for filenames (label.mop, label.out, ...) 

66 

67 Examples 

68 -------- 

69 

70 Use default values to do a single SCF calculation and print 

71 the forces (task='1SCF GRADIENTS'): 

72 

73 >>> from ase.build import molecule 

74 >>> from ase.calculators.mopac import MOPAC 

75 >>> 

76 >>> atoms = molecule('O2') 

77 >>> atoms.calc = MOPAC(label='O2') 

78 >>> atoms.get_potential_energy() 

79 >>> eigs = atoms.calc.get_eigenvalues() 

80 >>> somos = atoms.calc.get_somo_levels() 

81 >>> homo, lumo = atoms.calc.get_homo_lumo_levels() 

82 

83 Use the internal geometry optimization of Mopac: 

84 

85 >>> atoms = molecule('H2') 

86 >>> atoms.calc = MOPAC(label='H2', task='GRADIENTS') 

87 >>> atoms.get_potential_energy() 

88 

89 Read in and start from output file: 

90 

91 >>> atoms = MOPAC.read_atoms('H2') 

92 >>> atoms.calc.get_homo_lumo_levels() 

93 

94 """ 

95 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

96 label, atoms, **kwargs) 

97 

98 def write_input(self, atoms, properties=None, system_changes=None): 

99 FileIOCalculator.write_input(self, atoms, properties, system_changes) 

100 p = Parameters(self.parameters.copy()) 

101 

102 # Ensure DISP so total energy is available 

103 if 'DISP' not in p.task.split(): 

104 p.task = p.task + ' DISP' 

105 

106 # Build string to hold .mop input file: 

107 s = f'{p.method} {p.task} ' 

108 

109 if p.relscf: 

110 s += f'RELSCF={p.relscf} ' 

111 

112 # Write charge: 

113 if p.charge is None: 

114 charge = atoms.get_initial_charges().sum() 

115 else: 

116 charge = p.charge 

117 

118 if charge != 0: 

119 s += f'CHARGE={int(round(charge))} ' 

120 

121 magmom = int(round(abs(atoms.get_initial_magnetic_moments().sum()))) 

122 if magmom: 

123 s += (['DOUBLET', 'TRIPLET', 'QUARTET', 'QUINTET'][magmom - 1] + 

124 ' UHF ') 

125 

126 s += '\nTitle: ASE calculation\n\n' 

127 

128 # Write coordinates: 

129 for xyz, symbol in zip(atoms.positions, atoms.get_chemical_symbols()): 

130 s += ' {:2} {} 1 {} 1 {} 1\n'.format(symbol, *xyz) 

131 

132 for v, pbc in zip(atoms.cell, atoms.pbc): 

133 if pbc: 

134 s += 'Tv {} {} {}\n'.format(*v) 

135 

136 with open(self.label + '.mop', 'w') as fd: 

137 fd.write(s) 

138 

139 def get_spin_polarized(self): 

140 return self.nspins == 2 

141 

142 def get_index(self, lines, pattern): 

143 for i, line in enumerate(lines): 

144 if line.find(pattern) != -1: 

145 return i 

146 

147 def read(self, label): 

148 FileIOCalculator.read(self, label) 

149 if not os.path.isfile(self.label + '.out'): 

150 raise ReadError 

151 

152 with open(self.label + '.out') as fd: 

153 lines = fd.readlines() 

154 

155 self.parameters = Parameters(task='', method='') 

156 p = self.parameters 

157 parm_line = self.read_parameters_from_file(lines) 

158 for keyword in parm_line.split(): 

159 if 'RELSCF' in keyword: 

160 p.relscf = float(keyword.split('=')[-1]) 

161 elif keyword in self.methods: 

162 p.method = keyword 

163 else: 

164 p.task += keyword + ' ' 

165 

166 p.task = p.task.rstrip() 

167 if 'charge' not in p: 

168 p.charge = None 

169 

170 self.atoms = self.read_atoms_from_file(lines) 

171 self.read_results() 

172 

173 def read_atoms_from_file(self, lines): 

174 """Read the Atoms from the output file stored as list of str in lines. 

175 Parameters 

176 ---------- 

177 

178 lines: list of str 

179 """ 

180 # first try to read from final point (last image) 

181 i = self.get_index(lines, 'FINAL POINT AND DERIVATIVES') 

182 if i is None: # XXX should we read it from the input file? 

183 assert 0, 'Not implemented' 

184 

185 lines1 = lines[i:] 

186 i = self.get_index(lines1, 'CARTESIAN COORDINATES') 

187 j = i + 2 

188 symbols = [] 

189 positions = [] 

190 while not lines1[j].isspace(): # continue until we hit a blank line 

191 l = lines1[j].split() 

192 symbols.append(l[1]) 

193 positions.append([float(c) for c in l[2: 2 + 3]]) 

194 j += 1 

195 

196 return Atoms(symbols=symbols, positions=positions) 

197 

198 def read_parameters_from_file(self, lines): 

199 """Find and return the line that defines a Mopac calculation 

200 

201 Parameters 

202 ---------- 

203 

204 lines: list of str 

205 """ 

206 for i, line in enumerate(lines): 

207 if line.find('CALCULATION DONE:') != -1: 

208 break 

209 

210 lines1 = lines[i:] 

211 for i, line in enumerate(lines1): 

212 if line.find('****') != -1: 

213 return lines1[i + 1] 

214 

215 def read_results(self): 

216 """Read the results, such as energy, forces, eigenvalues, etc. 

217 """ 

218 FileIOCalculator.read(self, self.label) 

219 if not os.path.isfile(self.label + '.out'): 

220 raise ReadError 

221 

222 with open(self.label + '.out') as fd: 

223 lines = fd.readlines() 

224 

225 self.results['version'] = get_version_number(lines) 

226 

227 total_energy_regex = re.compile( 

228 r'^\s+TOTAL ENERGY\s+\=\s+(-?\d+\.\d+) EV') 

229 final_heat_regex = re.compile( 

230 r'^\s+FINAL HEAT OF FORMATION\s+\=\s+(-?\d+\.\d+) KCAL/MOL') 

231 

232 for i, line in enumerate(lines): 

233 if total_energy_regex.match(line): 

234 self.results['total_energy'] = float( 

235 total_energy_regex.match(line).groups()[0]) 

236 elif final_heat_regex.match(line): 

237 self.results['final_hof'] = float(final_heat_regex.match(line) 

238 .groups()[0]) * kcal / mol 

239 elif line.find('NO. OF FILLED LEVELS') != -1: 

240 self.nspins = 1 

241 self.no_occ_levels = int(line.split()[-1]) 

242 elif line.find('NO. OF ALPHA ELECTRON') != -1: 

243 self.nspins = 2 

244 self.no_alpha_electrons = int(line.split()[-1]) 

245 self.no_beta_electrons = int(lines[i + 1].split()[-1]) 

246 self.results['magmom'] = abs(self.no_alpha_electrons - 

247 self.no_beta_electrons) 

248 elif line.find('FINAL POINT AND DERIVATIVES') != -1: 

249 forces = [-float(line.split()[6]) 

250 for line in lines[i + 3:i + 3 + 3 * len(self.atoms)]] 

251 self.results['forces'] = np.array( 

252 forces).reshape((-1, 3)) * kcal / mol 

253 elif line.find('EIGENVALUES') != -1: 

254 if line.find('ALPHA') != -1: 

255 j = i + 1 

256 eigs_alpha = [] 

257 while not lines[j].isspace(): 

258 eigs_alpha += [float(eps) for eps in lines[j].split()] 

259 j += 1 

260 elif line.find('BETA') != -1: 

261 j = i + 1 

262 eigs_beta = [] 

263 while not lines[j].isspace(): 

264 eigs_beta += [float(eps) for eps in lines[j].split()] 

265 j += 1 

266 eigs = np.array([eigs_alpha, eigs_beta]).reshape(2, 1, -1) 

267 self.eigenvalues = eigs 

268 else: 

269 eigs = [] 

270 j = i + 1 

271 while not lines[j].isspace(): 

272 eigs += [float(e) for e in lines[j].split()] 

273 j += 1 

274 self.eigenvalues = np.array(eigs).reshape(1, 1, -1) 

275 elif line.find('DIPOLE ') != -1: 

276 self.results['dipole'] = np.array( 

277 lines[i + 3].split()[1:1 + 3], float) * Debye 

278 

279 # Developers recommend using final HOF as it includes dispersion etc. 

280 # For backward compatibility we won't change the results of old MOPAC 

281 # calculations... yet 

282 if version.parse(self.results['version']) >= version.parse('22'): 

283 self.results['energy'] = self.results['final_hof'] 

284 else: 

285 warn("Using a version of MOPAC lower than v22: ASE will use " 

286 "TOTAL ENERGY as the potential energy. In future, " 

287 "FINAL HEAT OF FORMATION will be preferred for all versions.") 

288 self.results['energy'] = self.results['total_energy'] 

289 self.results['free_energy'] = self.results['energy'] 

290 

291 def get_eigenvalues(self, kpt=0, spin=0): 

292 return self.eigenvalues[spin, kpt] 

293 

294 def get_homo_lumo_levels(self): 

295 eigs = self.eigenvalues 

296 if self.nspins == 1: 

297 nocc = self.no_occ_levels 

298 return np.array([eigs[0, 0, nocc - 1], eigs[0, 0, nocc]]) 

299 else: 

300 na = self.no_alpha_electrons 

301 nb = self.no_beta_electrons 

302 if na == 0: 

303 return None, self.eigenvalues[1, 0, nb - 1] 

304 elif nb == 0: 

305 return self.eigenvalues[0, 0, na - 1], None 

306 else: 

307 eah, eal = eigs[0, 0, na - 1: na + 1] 

308 ebh, ebl = eigs[1, 0, nb - 1: nb + 1] 

309 return np.array([max(eah, ebh), min(eal, ebl)]) 

310 

311 def get_somo_levels(self): 

312 assert self.nspins == 2 

313 na, nb = self.no_alpha_electrons, self.no_beta_electrons 

314 if na == 0: 

315 return None, self.eigenvalues[1, 0, nb - 1] 

316 elif nb == 0: 

317 return self.eigenvalues[0, 0, na - 1], None 

318 else: 

319 return np.array([self.eigenvalues[0, 0, na - 1], 

320 self.eigenvalues[1, 0, nb - 1]]) 

321 

322 def get_final_heat_of_formation(self): 

323 """Final heat of formation as reported in the Mopac output file""" 

324 warn("This method is deprecated, please use " 

325 "MOPAC.results['final_hof']", DeprecationWarning) 

326 return self.results['final_hof'] 

327 

328 @property 

329 def final_hof(self): 

330 warn("This property is deprecated, please use " 

331 "MOPAC.results['final_hof']", DeprecationWarning) 

332 return self.results['final_hof'] 

333 

334 @final_hof.setter 

335 def final_hof(self, new_hof): 

336 warn("This property is deprecated, please use " 

337 "MOPAC.results['final_hof']", DeprecationWarning) 

338 self.results['final_hof'] = new_hof