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 

15from packaging import version 

16import numpy as np 

17 

18from ase import Atoms 

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

20from ase.units import kcal, mol, Debye 

21 

22 

23class MOPAC(FileIOCalculator): 

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

25 'magmom', 'free_energy'] 

26 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 def __init__(self, restart=None, 

40 ignore_bad_restart_file=FileIOCalculator._deprecated, 

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

42 """Construct MOPAC-calculator object. 

43 

44 Parameters: 

45 

46 label: str 

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

48 

49 Examples: 

50 

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

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

53 

54 >>> from ase.build import molecule 

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

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

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

58 >>> atoms.get_potential_energy() 

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

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

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

62 

63 Use the internal geometry optimization of Mopac: 

64 

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

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

67 >>> atoms.get_potential_energy() 

68 

69 Read in and start from output file: 

70 

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

72 >>> atoms.calc.get_homo_lumo_levels() 

73 

74 """ 

75 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

76 label, atoms, **kwargs) 

77 

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

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

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

81 

82 # Ensure DISP so total energy is available 

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

84 p.task = p.task + ' DISP' 

85 

86 # Build string to hold .mop input file: 

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

88 

89 if p.relscf: 

90 s += 'RELSCF={0} '.format(p.relscf) 

91 

92 # Write charge: 

93 if p.charge is None: 

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

95 else: 

96 charge = p.charge 

97 

98 if charge != 0: 

99 s += 'CHARGE={0} '.format(int(round(charge))) 

100 

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

102 if magmom: 

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

104 ' UHF ') 

105 

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

107 

108 # Write coordinates: 

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

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

111 

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

113 if pbc: 

114 s += 'Tv {0} {1} {2}\n'.format(*v) 

115 

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

117 fd.write(s) 

118 

119 def get_spin_polarized(self): 

120 return self.nspins == 2 

121 

122 def get_index(self, lines, pattern): 

123 for i, line in enumerate(lines): 

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

125 return i 

126 

127 def read(self, label): 

128 FileIOCalculator.read(self, label) 

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

130 raise ReadError 

131 

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

133 lines = fd.readlines() 

134 

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

136 p = self.parameters 

137 parm_line = self.read_parameters_from_file(lines) 

138 for keyword in parm_line.split(): 

139 if 'RELSCF' in keyword: 

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

141 elif keyword in self.methods: 

142 p.method = keyword 

143 else: 

144 p.task += keyword + ' ' 

145 

146 p.task = p.task.rstrip() 

147 if 'charge' not in p: 

148 p.charge = None 

149 

150 self.atoms = self.read_atoms_from_file(lines) 

151 self.read_results() 

152 

153 def read_atoms_from_file(self, lines): 

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

155 Parameters: 

156 

157 lines: list of str 

158 """ 

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

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

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

162 assert 0, 'Not implemented' 

163 

164 lines1 = lines[i:] 

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

166 j = i + 2 

167 symbols = [] 

168 positions = [] 

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

170 l = lines1[j].split() 

171 symbols.append(l[1]) 

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

173 j += 1 

174 

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

176 

177 def read_parameters_from_file(self, lines): 

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

179 

180 Parameters: 

181 

182 lines: list of str 

183 """ 

184 for i, line in enumerate(lines): 

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

186 break 

187 

188 lines1 = lines[i:] 

189 for i, line in enumerate(lines1): 

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

191 return lines1[i + 1] 

192 

193 def read_results(self): 

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

195 """ 

196 FileIOCalculator.read(self, self.label) 

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

198 raise ReadError 

199 

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

201 lines = fd.readlines() 

202 

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

204 

205 total_energy_regex = re.compile( 

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

207 final_heat_regex = re.compile( 

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

209 

210 for i, line in enumerate(lines): 

211 if total_energy_regex.match(line): 

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

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

214 elif final_heat_regex.match(line): 

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

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

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

218 self.nspins = 1 

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

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

221 self.nspins = 2 

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

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

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

225 self.no_beta_electrons) 

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

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

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

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

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

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

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

233 j = i + 1 

234 eigs_alpha = [] 

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

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

237 j += 1 

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

239 j = i + 1 

240 eigs_beta = [] 

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

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

243 j += 1 

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

245 self.eigenvalues = eigs 

246 else: 

247 eigs = [] 

248 j = i + 1 

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

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

251 j += 1 

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

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

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

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

256 

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

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

259 # calculations... yet 

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

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

262 else: 

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

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

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

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

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

268 

269 @staticmethod 

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

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

272 for line in lines: 

273 match = version_regex.match(line) 

274 if match: 

275 return match.groups()[0] 

276 else: 

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

278 

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

280 return self.eigenvalues[spin, kpt] 

281 

282 def get_homo_lumo_levels(self): 

283 eigs = self.eigenvalues 

284 if self.nspins == 1: 

285 nocc = self.no_occ_levels 

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

287 else: 

288 na = self.no_alpha_electrons 

289 nb = self.no_beta_electrons 

290 if na == 0: 

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

292 elif nb == 0: 

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

294 else: 

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

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

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

298 

299 def get_somo_levels(self): 

300 assert self.nspins == 2 

301 na, nb = self.no_alpha_electrons, 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 return np.array([self.eigenvalues[0, 0, na - 1], 

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

309 

310 def get_final_heat_of_formation(self): 

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

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

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

314 return self.results['final_hof'] 

315 

316 @property 

317 def final_hof(self): 

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

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

320 return self.results['final_hof'] 

321 

322 @final_hof.setter 

323 def final_hof(self, new_hof): 

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

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

326 self.results['final_hof'] = new_hof