Coverage for ase / io / gamess_us.py: 75.39%

1"""Module for GAMESS US IO.""" 

2 

3from __future__ import annotations 

4 

5import os 

6import re 

7from copy import deepcopy 

8from io import TextIOBase 

9from subprocess import TimeoutExpired, call 

10 

11import numpy as np 

12 

13from ase import Atoms 

14from ase.calculators.singlepoint import SinglePointCalculator 

15from ase.units import Bohr, Debye, Hartree 

16from ase.utils import reader, workdir, writer 

17 

18 

19def _format_value(val: bool | str) -> str: 

20 if isinstance(val, bool): 

21 return '.t.' if val else '.f.' 

22 return str(val).upper() 

23 

24 

25def _write_block(name: str, args: dict[str, bool | str]) -> str: 

26 out = [f' ${name.upper()}'] 

27 for key, val in args.items(): 

28 out.append(f' {key.upper()}={_format_value(val)}') 

29 out.append(' $END') 

30 return '\n'.join(out) 

31 

32 

33def _write_geom(atoms: Atoms, basis_spec: dict[str | int, str] | None) -> str: 

34 out = [' $DATA', atoms.get_chemical_formula(), 'C1'] 

35 for i, atom in enumerate(atoms): 

36 out.append( 

37 '{:<3} {:>3} {:20.13e} {:20.13e} {:20.13e}'.format( 

38 atom.symbol, atom.number, *atom.position 

39 ) 

40 ) 

41 if basis_spec is not None: 

42 basis = basis_spec.get(i) 

43 if basis is None: 

44 basis = basis_spec.get(atom.symbol) 

45 if basis is None: 

46 msg = f'Could not find an appropriate basis set for atom {i}' 

47 raise ValueError(msg) 

48 out += [basis, ''] 

49 out.append(' $END') 

50 return '\n'.join(out) 

51 

52 

53def _write_ecp(atoms: Atoms, ecp: dict[int | str, str]) -> str: 

54 out = [' $ECP'] 

55 for i, symbol in enumerate(atoms.symbols): 

56 if i in ecp: 

57 out.append(ecp[i]) 

58 elif symbol in ecp: 

59 out.append(ecp[symbol]) 

60 else: 

61 msg = f'Could not find an appropriate ECP for atom {i}' 

62 raise ValueError(msg) 

63 out.append(' $END') 

64 return '\n'.join(out) 

65 

66 

67_xc = dict(LDA='SVWN') 

68 

69 

70@writer 

71def write_gamess_us_in(fd: TextIOBase, atoms: Atoms, properties=None, **params): 

72 params = deepcopy(params) 

73 

74 if properties is None: 

75 properties = ['energy'] 

76 

77 # set RUNTYP from properties iff value not provided by the user 

78 contrl = params.pop('contrl', {}) 

79 if 'runtyp' not in contrl: 

80 if 'forces' in properties: 

81 contrl['runtyp'] = 'gradient' 

82 else: 

83 contrl['runtyp'] = 'energy' 

84 

85 # Support xc keyword for functional specification 

86 xc = params.pop('xc', None) 

87 if xc is not None and 'dfttyp' not in contrl: 

88 contrl['dfttyp'] = _xc.get(xc.upper(), xc.upper()) 

89 

90 # Automatically determine multiplicity from magnetic moment 

91 magmom_tot = int(round(atoms.get_initial_magnetic_moments().sum())) 

92 if 'mult' not in contrl: 

93 contrl['mult'] = abs(magmom_tot) + 1 

94 

95 # Since we're automatically determining multiplicity, we also 

96 # need to automatically switch to UHF when the multiplicity 

97 # is not 1 

98 if 'scftyp' not in contrl: 

99 contrl['scftyp'] = 'rhf' if contrl['mult'] == 1 else 'uhf' 

100 

101 # effective core potentials 

102 ecp = params.pop('ecp', None) 

103 if ecp is not None and 'pp' not in contrl: 

104 contrl['pp'] = 'READ' 

105 

106 # If no basis set is provided, use 3-21G by default. 

107 basis_spec = None 

108 if 'basis' not in params: 

109 params['basis'] = dict(gbasis='N21', ngauss=3) 

110 else: 

111 keys = set(params['basis']) 

112 # Check if the user is specifying a literal per-atom basis. 

113 # We assume they are passing a per-atom basis if the keys of the 

114 # basis dict are atom symbols, or if they are atom indices, or 

115 # a mixture of both. 

116 if keys.intersection(set(atoms.symbols)) or any( 

117 map(lambda x: isinstance(x, int), keys) 

118 ): 

119 basis_spec = params.pop('basis') 

120 

121 out = [_write_block('contrl', contrl)] 

122 out += [_write_block(*item) for item in params.items()] 

123 out.append(_write_geom(atoms, basis_spec)) 

124 if ecp is not None: 

125 out.append(_write_ecp(atoms, ecp)) 

126 fd.write('\n\n'.join(out)) 

127 

128 

129_geom_re = re.compile(r'^\s*ATOM\s+ATOMIC\s+COORDINATES') 

130_atom_re = re.compile(r'^\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s*\n') 

131_energy_re = re.compile(r'^\s*FINAL [\S\s]+ ENERGY IS\s+(\S+) AFTER') 

132_grad_re = re.compile(r'^\s*GRADIENT OF THE ENERGY\s*') 

133_charges_re = re.compile(r'^\s*TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS\s*') 

134_dipole_re = re.compile(r'^\s+DX\s+DY\s+DZ\s+\/D\/\s+\(DEBYE\)') 

135 

136 

137@reader 

138def read_gamess_us_out(fd): 

139 atoms = None 

140 energy = None 

141 forces = None 

142 charges = None 

143 dipole = None 

144 for line in fd: 

145 # Geometry 

146 if _geom_re.match(line): 

147 fd.readline() 

148 symbols = [] 

149 pos = [] 

150 while True: 

151 atom = _atom_re.match(fd.readline()) 

152 if atom is None: 

153 break 

154 symbol, _, x, y, z = atom.groups() 

155 symbols.append(symbol.capitalize()) 

156 pos.append(list(map(float, [x, y, z]))) 

157 atoms = Atoms(symbols, np.array(pos) * Bohr) 

158 continue 

159 

160 # Energy 

161 ematch = _energy_re.match(line) 

162 if ematch is not None: 

163 energy = float(ematch.group(1)) * Hartree 

164 

165 # MPn energy. Supplants energy parsed above. 

166 elif line.strip().startswith('TOTAL ENERGY'): 

167 energy = float(line.strip().split()[-1]) * Hartree 

168 

169 # Higher-level energy (e.g. coupled cluster) 

170 # Supplants energies parsed above. 

171 elif line.strip().startswith('THE FOLLOWING METHOD AND ENERGY'): 

172 energy = float(fd.readline().strip().split()[-1]) * Hartree 

173 

174 elif _charges_re.match(line): 

175 fd.readline() 

176 charges = np.array( 

177 [float(fd.readline().split()[3]) for _ in symbols], 

178 ) 

179 

180 # Gradients 

181 elif _grad_re.match(line): 

182 for _ in range(3): 

183 fd.readline() 

184 grad = [] 

185 while True: 

186 atom = _atom_re.match(fd.readline()) 

187 if atom is None: 

188 break 

189 grad.append(list(map(float, atom.groups()[2:]))) 

190 forces = -np.array(grad) * Hartree / Bohr 

191 elif _dipole_re.match(line): 

192 dipole = np.array(list(map(float, fd.readline().split()[:3]))) 

193 dipole *= Debye 

194 

195 atoms.calc = SinglePointCalculator( 

196 atoms, 

197 energy=energy, 

198 forces=forces, 

199 charges=charges, 

200 dipole=dipole, 

201 ) 

202 return atoms 

203 

204 

205@reader 

206def read_gamess_us_punch(fd): 

207 atoms = None 

208 energy = None 

209 forces = None 

210 charges = None 

211 dipole = None 

212 for line in fd: 

213 if line.strip() == '$DATA': 

214 symbols = [] 

215 pos = [] 

216 while line.strip() != '$END': 

217 line = fd.readline() 

218 atom = _atom_re.match(line) 

219 if atom is None: 

220 # The basis set specification is interlaced with the 

221 # molecular geometry. We don't care about the basis 

222 # set, so ignore lines that don't match the pattern. 

223 continue 

224 symbols.append(atom.group(1).capitalize()) 

225 pos.append(list(map(float, atom.group(3, 4, 5)))) 

226 atoms = Atoms(symbols, np.array(pos)) 

227 elif line.startswith('E('): 

228 energy = float(line.split()[1][:-1]) * Hartree 

229 elif line.strip().startswith('POPULATION ANALYSIS'): 

230 # Mulliken charges 

231 charges = np.array( 

232 [float(fd.readline().split()[2]) for _ in symbols], 

233 ) 

234 elif line.strip().startswith('DIPOLE'): 

235 dipole = np.array(list(map(float, line.split()[1:]))) * Debye 

236 elif line.strip() == '$GRAD': 

237 # The gradient block also contains the energy, which we prefer 

238 # over the energy obtained above because it is more likely to 

239 # be consistent with the gradients. It probably doesn't actually 

240 # make a difference though. 

241 energy = float(fd.readline().split()[1]) * Hartree 

242 grad = [] 

243 while line.strip() != '$END': 

244 line = fd.readline() 

245 atom = _atom_re.match(line) 

246 if atom is None: 

247 continue 

248 grad.append(list(map(float, atom.group(3, 4, 5)))) 

249 forces = -np.array(grad) * Hartree / Bohr 

250 

251 atoms.calc = SinglePointCalculator( 

252 atoms, 

253 energy=energy, 

254 forces=forces, 

255 charges=charges, 

256 dipole=dipole, 

257 ) 

258 

259 return atoms 

260 

261 

262def clean_userscr(userscr: str, prefix: str) -> None: 

263 for fname in os.listdir(userscr): 

264 tokens = fname.split('.') 

265 if tokens[0] == prefix and tokens[-1] != 'bak': 

266 fold = os.path.join(userscr, fname) 

267 os.rename(fold, fold + '.bak') 

268 

269 

270def get_userscr(prefix: str, command: str) -> str | None: 

271 prefix_test = prefix + '_test' 

272 command = command.replace('PREFIX', prefix_test) 

273 with workdir(prefix_test, mkdir=True): 

274 try: 

275 call(command, shell=True, timeout=2) 

276 except TimeoutExpired: 

277 pass 

278 

279 try: 

280 with open(f'{prefix_test}.log', encoding='utf-8') as fd: 

281 for line in fd: 

282 if line.startswith('GAMESS supplementary output files'): 

283 return ' '.join(line.split(' ')[8:]).strip() 

284 except FileNotFoundError: 

285 return None 

286 

287 return None