1import os 

2import re 

3from subprocess import call, TimeoutExpired 

4from copy import deepcopy 

5 

6import numpy as np 

7 

8from ase import Atoms 

9from ase.utils import workdir 

10from ase.units import Hartree, Bohr, Debye 

11from ase.calculators.singlepoint import SinglePointCalculator 

12 

13 

14def _format_value(val): 

15 if isinstance(val, bool): 

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

17 return str(val).upper() 

18 

19 

20def _write_block(name, args): 

21 out = [' ${}'.format(name.upper())] 

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

23 out.append(' {}={}'.format(key.upper(), _format_value(val))) 

24 out.append(' $END') 

25 return '\n'.join(out) 

26 

27 

28def _write_geom(atoms, basis_spec): 

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

30 for i, atom in enumerate(atoms): 

31 out.append('{:<3} {:>3} {:20.13e} {:20.13e} {:20.13e}' 

32 .format(atom.symbol, atom.number, *atom.position)) 

33 if basis_spec is not None: 

34 basis = basis_spec.get(i) 

35 if basis is None: 

36 basis = basis_spec.get(atom.symbol) 

37 if basis is None: 

38 raise ValueError('Could not find an appropriate basis set ' 

39 'for atom number {}!'.format(i)) 

40 out += [basis, ''] 

41 out.append(' $END') 

42 return '\n'.join(out) 

43 

44 

45def _write_ecp(atoms, ecp): 

46 out = [' $ECP'] 

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

48 if i in ecp: 

49 out.append(ecp[i]) 

50 elif symbol in ecp: 

51 out.append(ecp[symbol]) 

52 else: 

53 raise ValueError('Could not find an appropriate ECP for ' 

54 'atom number {}!'.format(i)) 

55 out.append(' $END') 

56 return '\n'.join(out) 

57 

58 

59_xc = dict(LDA='SVWN') 

60 

61 

62def write_gamess_us_in(fd, atoms, properties=None, **params): 

63 params = deepcopy(params) 

64 

65 if properties is None: 

66 properties = ['energy'] 

67 

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

69 contrl = params.pop('contrl', dict()) 

70 if 'runtyp' not in contrl: 

71 if 'forces' in properties: 

72 contrl['runtyp'] = 'gradient' 

73 else: 

74 contrl['runtyp'] = 'energy' 

75 

76 # Support xc keyword for functional specification 

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

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

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

80 

81 # Automatically determine multiplicity from magnetic moment 

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

83 if 'mult' not in contrl: 

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

85 

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

87 # need to automatically switch to UHF when the multiplicity 

88 # is not 1 

89 if 'scftyp' not in contrl: 

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

91 

92 # effective core potentials 

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

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

95 contrl['pp'] = 'READ' 

96 

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

98 basis_spec = None 

99 if 'basis' not in params: 

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

101 else: 

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

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

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

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

106 # a mixture of both. 

107 if (keys.intersection(set(atoms.symbols)) 

108 or any(map(lambda x: isinstance(x, int), keys))): 

109 basis_spec = params.pop('basis') 

110 

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

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

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

114 if ecp is not None: 

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

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

117 

118 

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

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

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

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

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

124 

125 

126def read_gamess_us_out(fd): 

127 atoms = None 

128 energy = None 

129 forces = None 

130 dipole = None 

131 for line in fd: 

132 # Geometry 

133 if _geom_re.match(line): 

134 fd.readline() 

135 symbols = [] 

136 pos = [] 

137 while True: 

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

139 if atom is None: 

140 break 

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

142 symbols.append(symbol.capitalize()) 

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

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

145 continue 

146 

147 # Energy 

148 ematch = _energy_re.match(line) 

149 if ematch is not None: 

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

151 

152 # MPn energy. Supplants energy parsed above. 

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

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

155 

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

157 # Supplants energies parsed above. 

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

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

160 

161 # Gradients 

162 elif _grad_re.match(line): 

163 for _ in range(3): 

164 fd.readline() 

165 grad = [] 

166 while True: 

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

168 if atom is None: 

169 break 

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

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

172 elif _dipole_re.match(line): 

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

174 dipole *= Debye 

175 

176 atoms.calc = SinglePointCalculator(atoms, energy=energy, 

177 forces=forces, dipole=dipole) 

178 return atoms 

179 

180 

181def read_gamess_us_punch(fd): 

182 atoms = None 

183 energy = None 

184 forces = None 

185 dipole = None 

186 for line in fd: 

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

188 symbols = [] 

189 pos = [] 

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

191 line = fd.readline() 

192 atom = _atom_re.match(line) 

193 if atom is None: 

194 # The basis set specification is interlaced with the 

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

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

197 continue 

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

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

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

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

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

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

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

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

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

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

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

209 # make a difference though. 

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

211 grad = [] 

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

213 line = fd.readline() 

214 atom = _atom_re.match(line) 

215 if atom is None: 

216 continue 

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

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

219 

220 atoms.calc = SinglePointCalculator(atoms, energy=energy, forces=forces, 

221 dipole=dipole) 

222 

223 return atoms 

224 

225 

226def clean_userscr(userscr, prefix): 

227 for fname in os.listdir(userscr): 

228 tokens = fname.split('.') 

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

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

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

232 

233 

234def get_userscr(prefix, command): 

235 prefix_test = prefix + '_test' 

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

237 with workdir(prefix_test, mkdir=True): 

238 try: 

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

240 except TimeoutExpired: 

241 pass 

242 

243 try: 

244 with open(prefix_test + '.log') as fd: 

245 for line in fd: 

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

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

248 except FileNotFoundError: 

249 return None 

250 

251 return None