Coverage for /builds/kinetik161/ase/ase/calculators/mopac.py: 85.48%

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 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 >>> 

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

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

59 >>> atoms.get_potential_energy() 

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

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

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

63 

64 Use the internal geometry optimization of Mopac: 

65 

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

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

68 >>> atoms.get_potential_energy() 

69 

70 Read in and start from output file: 

71 

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

73 >>> atoms.calc.get_homo_lumo_levels() 

74 

75 """ 

76 FileIOCalculator.__init__(self, restart, ignore_bad_restart_file, 

77 label, atoms, **kwargs) 

78 

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

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

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

82 

83 # Ensure DISP so total energy is available 

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

85 p.task = p.task + ' DISP' 

86 

87 # Build string to hold .mop input file: 

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

89 

90 if p.relscf: 

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

92 

93 # Write charge: 

94 if p.charge is None: 

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

96 else: 

97 charge = p.charge 

98 

99 if charge != 0: 

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

101 

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

103 if magmom: 

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

105 ' UHF ') 

106 

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

108 

109 # Write coordinates: 

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

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

112 

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

114 if pbc: 

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

116 

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

118 fd.write(s) 

119 

120 def get_spin_polarized(self): 

121 return self.nspins == 2 

122 

123 def get_index(self, lines, pattern): 

124 for i, line in enumerate(lines): 

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

126 return i 

127 

128 def read(self, label): 

129 FileIOCalculator.read(self, label) 

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

131 raise ReadError 

132 

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

134 lines = fd.readlines() 

135 

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

137 p = self.parameters 

138 parm_line = self.read_parameters_from_file(lines) 

139 for keyword in parm_line.split(): 

140 if 'RELSCF' in keyword: 

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

142 elif keyword in self.methods: 

143 p.method = keyword 

144 else: 

145 p.task += keyword + ' ' 

146 

147 p.task = p.task.rstrip() 

148 if 'charge' not in p: 

149 p.charge = None 

150 

151 self.atoms = self.read_atoms_from_file(lines) 

152 self.read_results() 

153 

154 def read_atoms_from_file(self, lines): 

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

156 Parameters: 

157 

158 lines: list of str 

159 """ 

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

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

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

163 assert 0, 'Not implemented' 

164 

165 lines1 = lines[i:] 

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

167 j = i + 2 

168 symbols = [] 

169 positions = [] 

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

171 l = lines1[j].split() 

172 symbols.append(l[1]) 

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

174 j += 1 

175 

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

177 

178 def read_parameters_from_file(self, lines): 

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

180 

181 Parameters: 

182 

183 lines: list of str 

184 """ 

185 for i, line in enumerate(lines): 

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

187 break 

188 

189 lines1 = lines[i:] 

190 for i, line in enumerate(lines1): 

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

192 return lines1[i + 1] 

193 

194 def read_results(self): 

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

196 """ 

197 FileIOCalculator.read(self, self.label) 

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

199 raise ReadError 

200 

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

202 lines = fd.readlines() 

203 

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

205 

206 total_energy_regex = re.compile( 

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

208 final_heat_regex = re.compile( 

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

210 

211 for i, line in enumerate(lines): 

212 if total_energy_regex.match(line): 

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

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

215 elif final_heat_regex.match(line): 

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

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

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

219 self.nspins = 1 

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

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

222 self.nspins = 2 

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

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

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

226 self.no_beta_electrons) 

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

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

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

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

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

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

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

234 j = i + 1 

235 eigs_alpha = [] 

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

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

238 j += 1 

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

240 j = i + 1 

241 eigs_beta = [] 

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

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

244 j += 1 

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

246 self.eigenvalues = eigs 

247 else: 

248 eigs = [] 

249 j = i + 1 

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

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

252 j += 1 

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

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

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

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

257 

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

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

260 # calculations... yet 

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

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

263 else: 

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

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

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

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

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

269 

270 @staticmethod 

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

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

273 for line in lines: 

274 match = version_regex.match(line) 

275 if match: 

276 return match.groups()[0] 

277 else: 

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

279 

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

281 return self.eigenvalues[spin, kpt] 

282 

283 def get_homo_lumo_levels(self): 

284 eigs = self.eigenvalues 

285 if self.nspins == 1: 

286 nocc = self.no_occ_levels 

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

288 else: 

289 na = self.no_alpha_electrons 

290 nb = self.no_beta_electrons 

291 if na == 0: 

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

293 elif nb == 0: 

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

295 else: 

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

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

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

299 

300 def get_somo_levels(self): 

301 assert self.nspins == 2 

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

303 if na == 0: 

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

305 elif nb == 0: 

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

307 else: 

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

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

310 

311 def get_final_heat_of_formation(self): 

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

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

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

315 return self.results['final_hof'] 

316 

317 @property 

318 def final_hof(self): 

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

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

321 return self.results['final_hof'] 

322 

323 @final_hof.setter 

324 def final_hof(self, new_hof): 

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

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

327 self.results['final_hof'] = new_hof