Coverage for /builds/ase/ase/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 typing 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 label: str 

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

65 

66 Examples: 

67 

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

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

70 

71 >>> from ase.build import molecule 

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

73 >>> 

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

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

76 >>> atoms.get_potential_energy() 

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

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

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

80 

81 Use the internal geometry optimization of Mopac: 

82 

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

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

85 >>> atoms.get_potential_energy() 

86 

87 Read in and start from output file: 

88 

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

90 >>> atoms.calc.get_homo_lumo_levels() 

91 

92 """ 

93 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

94 label, atoms, **kwargs) 

95 

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

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

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

99 

100 # Ensure DISP so total energy is available 

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

102 p.task = p.task + ' DISP' 

103 

104 # Build string to hold .mop input file: 

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

106 

107 if p.relscf: 

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

109 

110 # Write charge: 

111 if p.charge is None: 

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

113 else: 

114 charge = p.charge 

115 

116 if charge != 0: 

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

118 

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

120 if magmom: 

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

122 ' UHF ') 

123 

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

125 

126 # Write coordinates: 

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

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

129 

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

131 if pbc: 

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

133 

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

135 fd.write(s) 

136 

137 def get_spin_polarized(self): 

138 return self.nspins == 2 

139 

140 def get_index(self, lines, pattern): 

141 for i, line in enumerate(lines): 

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

143 return i 

144 

145 def read(self, label): 

146 FileIOCalculator.read(self, label) 

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

148 raise ReadError 

149 

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

151 lines = fd.readlines() 

152 

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

154 p = self.parameters 

155 parm_line = self.read_parameters_from_file(lines) 

156 for keyword in parm_line.split(): 

157 if 'RELSCF' in keyword: 

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

159 elif keyword in self.methods: 

160 p.method = keyword 

161 else: 

162 p.task += keyword + ' ' 

163 

164 p.task = p.task.rstrip() 

165 if 'charge' not in p: 

166 p.charge = None 

167 

168 self.atoms = self.read_atoms_from_file(lines) 

169 self.read_results() 

170 

171 def read_atoms_from_file(self, lines): 

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

173 Parameters: 

174 

175 lines: list of str 

176 """ 

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

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

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

180 assert 0, 'Not implemented' 

181 

182 lines1 = lines[i:] 

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

184 j = i + 2 

185 symbols = [] 

186 positions = [] 

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

188 l = lines1[j].split() 

189 symbols.append(l[1]) 

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

191 j += 1 

192 

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

194 

195 def read_parameters_from_file(self, lines): 

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

197 

198 Parameters: 

199 

200 lines: list of str 

201 """ 

202 for i, line in enumerate(lines): 

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

204 break 

205 

206 lines1 = lines[i:] 

207 for i, line in enumerate(lines1): 

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

209 return lines1[i + 1] 

210 

211 def read_results(self): 

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

213 """ 

214 FileIOCalculator.read(self, self.label) 

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

216 raise ReadError 

217 

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

219 lines = fd.readlines() 

220 

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

222 

223 total_energy_regex = re.compile( 

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

225 final_heat_regex = re.compile( 

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

227 

228 for i, line in enumerate(lines): 

229 if total_energy_regex.match(line): 

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

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

232 elif final_heat_regex.match(line): 

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

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

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

236 self.nspins = 1 

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

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

239 self.nspins = 2 

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

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

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

243 self.no_beta_electrons) 

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

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

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

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

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

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

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

251 j = i + 1 

252 eigs_alpha = [] 

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

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

255 j += 1 

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

257 j = i + 1 

258 eigs_beta = [] 

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

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

261 j += 1 

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

263 self.eigenvalues = eigs 

264 else: 

265 eigs = [] 

266 j = i + 1 

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

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

269 j += 1 

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

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

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

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

274 

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

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

277 # calculations... yet 

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

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

280 else: 

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

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

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

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

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

286 

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

288 return self.eigenvalues[spin, kpt] 

289 

290 def get_homo_lumo_levels(self): 

291 eigs = self.eigenvalues 

292 if self.nspins == 1: 

293 nocc = self.no_occ_levels 

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

295 else: 

296 na = self.no_alpha_electrons 

297 nb = self.no_beta_electrons 

298 if na == 0: 

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

300 elif nb == 0: 

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

302 else: 

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

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

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

306 

307 def get_somo_levels(self): 

308 assert self.nspins == 2 

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

310 if na == 0: 

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

312 elif nb == 0: 

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

314 else: 

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

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

317 

318 def get_final_heat_of_formation(self): 

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

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

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

322 return self.results['final_hof'] 

323 

324 @property 

325 def final_hof(self): 

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

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

328 return self.results['final_hof'] 

329 

330 @final_hof.setter 

331 def final_hof(self, new_hof): 

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

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

334 self.results['final_hof'] = new_hof