Coverage for /builds/hweiske/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 

23class MOPAC(FileIOCalculator): 

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

25 'magmom', 'free_energy'] 

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

27 discard_results_on_any_change = True 

28 

29 default_parameters = dict( 

30 method='PM7', 

31 task='1SCF GRADIENTS', 

32 charge=None, 

33 relscf=0.0001) 

34 

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

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

37 'PM7-TS', 'RM1'] 

38 

39 fileio_rules = FileIOCalculator.ruleset( 

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

41 stdout_name='{prefix}.out') 

42 

43 def __init__(self, restart=None, 

44 ignore_bad_restart_file=FileIOCalculator._deprecated, 

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

46 """Construct MOPAC-calculator object. 

47 

48 Parameters: 

49 

50 label: str 

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

52 

53 Examples: 

54 

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

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

57 

58 >>> from ase.build import molecule 

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

60 >>> 

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

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

63 >>> atoms.get_potential_energy() 

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

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

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

67 

68 Use the internal geometry optimization of Mopac: 

69 

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

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

72 >>> atoms.get_potential_energy() 

73 

74 Read in and start from output file: 

75 

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

77 >>> atoms.calc.get_homo_lumo_levels() 

78 

79 """ 

80 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

81 label, atoms, **kwargs) 

82 

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

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

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

86 

87 # Ensure DISP so total energy is available 

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

89 p.task = p.task + ' DISP' 

90 

91 # Build string to hold .mop input file: 

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

93 

94 if p.relscf: 

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

96 

97 # Write charge: 

98 if p.charge is None: 

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

100 else: 

101 charge = p.charge 

102 

103 if charge != 0: 

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

105 

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

107 if magmom: 

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

109 ' UHF ') 

110 

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

112 

113 # Write coordinates: 

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

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

116 

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

118 if pbc: 

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

120 

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

122 fd.write(s) 

123 

124 def get_spin_polarized(self): 

125 return self.nspins == 2 

126 

127 def get_index(self, lines, pattern): 

128 for i, line in enumerate(lines): 

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

130 return i 

131 

132 def read(self, label): 

133 FileIOCalculator.read(self, label) 

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

135 raise ReadError 

136 

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

138 lines = fd.readlines() 

139 

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

141 p = self.parameters 

142 parm_line = self.read_parameters_from_file(lines) 

143 for keyword in parm_line.split(): 

144 if 'RELSCF' in keyword: 

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

146 elif keyword in self.methods: 

147 p.method = keyword 

148 else: 

149 p.task += keyword + ' ' 

150 

151 p.task = p.task.rstrip() 

152 if 'charge' not in p: 

153 p.charge = None 

154 

155 self.atoms = self.read_atoms_from_file(lines) 

156 self.read_results() 

157 

158 def read_atoms_from_file(self, lines): 

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

160 Parameters: 

161 

162 lines: list of str 

163 """ 

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

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

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

167 assert 0, 'Not implemented' 

168 

169 lines1 = lines[i:] 

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

171 j = i + 2 

172 symbols = [] 

173 positions = [] 

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

175 l = lines1[j].split() 

176 symbols.append(l[1]) 

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

178 j += 1 

179 

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

181 

182 def read_parameters_from_file(self, lines): 

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

184 

185 Parameters: 

186 

187 lines: list of str 

188 """ 

189 for i, line in enumerate(lines): 

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

191 break 

192 

193 lines1 = lines[i:] 

194 for i, line in enumerate(lines1): 

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

196 return lines1[i + 1] 

197 

198 def read_results(self): 

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

200 """ 

201 FileIOCalculator.read(self, self.label) 

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

203 raise ReadError 

204 

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

206 lines = fd.readlines() 

207 

208 self.results['version'] = self.get_version_from_file(lines) 

209 

210 total_energy_regex = re.compile( 

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

212 final_heat_regex = re.compile( 

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

214 

215 for i, line in enumerate(lines): 

216 if total_energy_regex.match(line): 

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

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

219 elif final_heat_regex.match(line): 

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

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

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

223 self.nspins = 1 

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

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

226 self.nspins = 2 

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

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

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

230 self.no_beta_electrons) 

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

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

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

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

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

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

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

238 j = i + 1 

239 eigs_alpha = [] 

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

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

242 j += 1 

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

244 j = i + 1 

245 eigs_beta = [] 

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

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

248 j += 1 

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

250 self.eigenvalues = eigs 

251 else: 

252 eigs = [] 

253 j = i + 1 

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

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

256 j += 1 

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

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

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

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

261 

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

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

264 # calculations... yet 

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

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

267 else: 

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

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

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

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

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

273 

274 @staticmethod 

275 def get_version_from_file(lines: Sequence[str]): 

276 version_regex = re.compile(r'^ \*\*\s+MOPAC (v[\.\d]+)\s+\*\*\s$') 

277 for line in lines: 

278 match = version_regex.match(line) 

279 if match: 

280 return match.groups()[0] 

281 return ValueError('Version number was not found in MOPAC output') 

282 

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

284 return self.eigenvalues[spin, kpt] 

285 

286 def get_homo_lumo_levels(self): 

287 eigs = self.eigenvalues 

288 if self.nspins == 1: 

289 nocc = self.no_occ_levels 

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

291 else: 

292 na = self.no_alpha_electrons 

293 nb = self.no_beta_electrons 

294 if na == 0: 

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

296 elif nb == 0: 

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

298 else: 

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

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

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

302 

303 def get_somo_levels(self): 

304 assert self.nspins == 2 

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

306 if na == 0: 

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

308 elif nb == 0: 

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

310 else: 

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

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

313 

314 def get_final_heat_of_formation(self): 

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

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

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

318 return self.results['final_hof'] 

319 

320 @property 

321 def final_hof(self): 

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

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

324 return self.results['final_hof'] 

325 

326 @final_hof.setter 

327 def final_hof(self, new_hof): 

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

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

330 self.results['final_hof'] = new_hof