Coverage for /builds/debichem-team/python-ase/ase/calculators/mopac.py: 85.56%

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

2 

3Set $ASE_MOPAC_COMMAND to something like:: 

4 

5 LD_LIBRARY_PATH=/path/to/lib/ \ 

6 MOPAC_LICENSE=/path/to/license \ 

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

8 

9""" 

10import os 

11import re 

12from typing import Sequence 

13from warnings import warn 

14 

15import numpy as np 

16from packaging import version 

17 

18from ase import Atoms 

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

20from ase.units import Debye, kcal, mol 

21 

22 

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

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

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

26 

27 for line in lines[:10]: 

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

29 if match: 

30 return match.group(1) 

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

32 

33 

34class MOPAC(FileIOCalculator): 

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

36 'magmom', 'free_energy'] 

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

38 discard_results_on_any_change = True 

39 

40 default_parameters = dict( 

41 method='PM7', 

42 task='1SCF GRADIENTS', 

43 charge=None, 

44 relscf=0.0001) 

45 

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

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

48 'PM7-TS', 'RM1'] 

49 

50 fileio_rules = FileIOCalculator.ruleset( 

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

52 stdout_name='{prefix}.out') 

53 

54 def __init__(self, restart=None, 

55 ignore_bad_restart_file=FileIOCalculator._deprecated, 

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

57 """Construct MOPAC-calculator object. 

58 

59 Parameters: 

60 

61 label: str 

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

63 

64 Examples: 

65 

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

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

68 

69 >>> from ase.build import molecule 

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

71 >>> 

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

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

74 >>> atoms.get_potential_energy() 

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

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

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

78 

79 Use the internal geometry optimization of Mopac: 

80 

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

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

83 >>> atoms.get_potential_energy() 

84 

85 Read in and start from output file: 

86 

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

88 >>> atoms.calc.get_homo_lumo_levels() 

89 

90 """ 

91 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

92 label, atoms, **kwargs) 

93 

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

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

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

97 

98 # Ensure DISP so total energy is available 

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

100 p.task = p.task + ' DISP' 

101 

102 # Build string to hold .mop input file: 

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

104 

105 if p.relscf: 

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

107 

108 # Write charge: 

109 if p.charge is None: 

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

111 else: 

112 charge = p.charge 

113 

114 if charge != 0: 

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

116 

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

118 if magmom: 

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

120 ' UHF ') 

121 

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

123 

124 # Write coordinates: 

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

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

127 

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

129 if pbc: 

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

131 

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

133 fd.write(s) 

134 

135 def get_spin_polarized(self): 

136 return self.nspins == 2 

137 

138 def get_index(self, lines, pattern): 

139 for i, line in enumerate(lines): 

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

141 return i 

142 

143 def read(self, label): 

144 FileIOCalculator.read(self, label) 

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

146 raise ReadError 

147 

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

149 lines = fd.readlines() 

150 

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

152 p = self.parameters 

153 parm_line = self.read_parameters_from_file(lines) 

154 for keyword in parm_line.split(): 

155 if 'RELSCF' in keyword: 

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

157 elif keyword in self.methods: 

158 p.method = keyword 

159 else: 

160 p.task += keyword + ' ' 

161 

162 p.task = p.task.rstrip() 

163 if 'charge' not in p: 

164 p.charge = None 

165 

166 self.atoms = self.read_atoms_from_file(lines) 

167 self.read_results() 

168 

169 def read_atoms_from_file(self, lines): 

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

171 Parameters: 

172 

173 lines: list of str 

174 """ 

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

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

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

178 assert 0, 'Not implemented' 

179 

180 lines1 = lines[i:] 

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

182 j = i + 2 

183 symbols = [] 

184 positions = [] 

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

186 l = lines1[j].split() 

187 symbols.append(l[1]) 

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

189 j += 1 

190 

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

192 

193 def read_parameters_from_file(self, lines): 

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

195 

196 Parameters: 

197 

198 lines: list of str 

199 """ 

200 for i, line in enumerate(lines): 

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

202 break 

203 

204 lines1 = lines[i:] 

205 for i, line in enumerate(lines1): 

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

207 return lines1[i + 1] 

208 

209 def read_results(self): 

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

211 """ 

212 FileIOCalculator.read(self, self.label) 

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

214 raise ReadError 

215 

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

217 lines = fd.readlines() 

218 

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

220 

221 total_energy_regex = re.compile( 

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

223 final_heat_regex = re.compile( 

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

225 

226 for i, line in enumerate(lines): 

227 if total_energy_regex.match(line): 

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

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

230 elif final_heat_regex.match(line): 

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

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

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

234 self.nspins = 1 

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

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

237 self.nspins = 2 

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

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

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

241 self.no_beta_electrons) 

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

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

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

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

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

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

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

249 j = i + 1 

250 eigs_alpha = [] 

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

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

253 j += 1 

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

255 j = i + 1 

256 eigs_beta = [] 

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

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

259 j += 1 

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

261 self.eigenvalues = eigs 

262 else: 

263 eigs = [] 

264 j = i + 1 

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

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

267 j += 1 

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

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

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

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

272 

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

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

275 # calculations... yet 

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

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

278 else: 

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

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

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

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

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

284 

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

286 return self.eigenvalues[spin, kpt] 

287 

288 def get_homo_lumo_levels(self): 

289 eigs = self.eigenvalues 

290 if self.nspins == 1: 

291 nocc = self.no_occ_levels 

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

293 else: 

294 na = self.no_alpha_electrons 

295 nb = self.no_beta_electrons 

296 if na == 0: 

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

298 elif nb == 0: 

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

300 else: 

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

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

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

304 

305 def get_somo_levels(self): 

306 assert self.nspins == 2 

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

308 if na == 0: 

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

310 elif nb == 0: 

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

312 else: 

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

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

315 

316 def get_final_heat_of_formation(self): 

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

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

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

320 return self.results['final_hof'] 

321 

322 @property 

323 def final_hof(self): 

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

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

326 return self.results['final_hof'] 

327 

328 @final_hof.setter 

329 def final_hof(self, new_hof): 

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

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

332 self.results['final_hof'] = new_hof