Coverage for /builds/ase/ase/ase/io/gamess_us.py: 34.66%
176 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import os
4import re
5from copy import deepcopy
6from subprocess import TimeoutExpired, call
8import numpy as np
10from ase import Atoms
11from ase.calculators.singlepoint import SinglePointCalculator
12from ase.units import Bohr, Debye, Hartree
13from ase.utils import workdir
16def _format_value(val):
17 if isinstance(val, bool):
18 return '.t.' if val else '.f.'
19 return str(val).upper()
22def _write_block(name, args):
23 out = [f' ${name.upper()}']
24 for key, val in args.items():
25 out.append(f' {key.upper()}={_format_value(val)}')
26 out.append(' $END')
27 return '\n'.join(out)
30def _write_geom(atoms, basis_spec):
31 out = [' $DATA', atoms.get_chemical_formula(), 'C1']
32 for i, atom in enumerate(atoms):
33 out.append('{:<3} {:>3} {:20.13e} {:20.13e} {:20.13e}'
34 .format(atom.symbol, atom.number, *atom.position))
35 if basis_spec is not None:
36 basis = basis_spec.get(i)
37 if basis is None:
38 basis = basis_spec.get(atom.symbol)
39 if basis is None:
40 raise ValueError('Could not find an appropriate basis set '
41 'for atom number {}!'.format(i))
42 out += [basis, '']
43 out.append(' $END')
44 return '\n'.join(out)
47def _write_ecp(atoms, ecp):
48 out = [' $ECP']
49 for i, symbol in enumerate(atoms.symbols):
50 if i in ecp:
51 out.append(ecp[i])
52 elif symbol in ecp:
53 out.append(ecp[symbol])
54 else:
55 raise ValueError('Could not find an appropriate ECP for '
56 'atom number {}!'.format(i))
57 out.append(' $END')
58 return '\n'.join(out)
61_xc = dict(LDA='SVWN')
64def write_gamess_us_in(fd, atoms, properties=None, **params):
65 params = deepcopy(params)
67 if properties is None:
68 properties = ['energy']
70 # set RUNTYP from properties iff value not provided by the user
71 contrl = params.pop('contrl', {})
72 if 'runtyp' not in contrl:
73 if 'forces' in properties:
74 contrl['runtyp'] = 'gradient'
75 else:
76 contrl['runtyp'] = 'energy'
78 # Support xc keyword for functional specification
79 xc = params.pop('xc', None)
80 if xc is not None and 'dfttyp' not in contrl:
81 contrl['dfttyp'] = _xc.get(xc.upper(), xc.upper())
83 # Automatically determine multiplicity from magnetic moment
84 magmom_tot = int(round(atoms.get_initial_magnetic_moments().sum()))
85 if 'mult' not in contrl:
86 contrl['mult'] = abs(magmom_tot) + 1
88 # Since we're automatically determining multiplicity, we also
89 # need to automatically switch to UHF when the multiplicity
90 # is not 1
91 if 'scftyp' not in contrl:
92 contrl['scftyp'] = 'rhf' if contrl['mult'] == 1 else 'uhf'
94 # effective core potentials
95 ecp = params.pop('ecp', None)
96 if ecp is not None and 'pp' not in contrl:
97 contrl['pp'] = 'READ'
99 # If no basis set is provided, use 3-21G by default.
100 basis_spec = None
101 if 'basis' not in params:
102 params['basis'] = dict(gbasis='N21', ngauss=3)
103 else:
104 keys = set(params['basis'])
105 # Check if the user is specifying a literal per-atom basis.
106 # We assume they are passing a per-atom basis if the keys of the
107 # basis dict are atom symbols, or if they are atom indices, or
108 # a mixture of both.
109 if (keys.intersection(set(atoms.symbols))
110 or any(map(lambda x: isinstance(x, int), keys))):
111 basis_spec = params.pop('basis')
113 out = [_write_block('contrl', contrl)]
114 out += [_write_block(*item) for item in params.items()]
115 out.append(_write_geom(atoms, basis_spec))
116 if ecp is not None:
117 out.append(_write_ecp(atoms, ecp))
118 fd.write('\n\n'.join(out))
121_geom_re = re.compile(r'^\s*ATOM\s+ATOMIC\s+COORDINATES')
122_atom_re = re.compile(r'^\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s*\n')
123_energy_re = re.compile(r'^\s*FINAL [\S\s]+ ENERGY IS\s+(\S+) AFTER')
124_grad_re = re.compile(r'^\s*GRADIENT OF THE ENERGY\s*')
125_dipole_re = re.compile(r'^\s+DX\s+DY\s+DZ\s+\/D\/\s+\(DEBYE\)')
128def read_gamess_us_out(fd):
129 atoms = None
130 energy = None
131 forces = None
132 dipole = None
133 for line in fd:
134 # Geometry
135 if _geom_re.match(line):
136 fd.readline()
137 symbols = []
138 pos = []
139 while True:
140 atom = _atom_re.match(fd.readline())
141 if atom is None:
142 break
143 symbol, _, x, y, z = atom.groups()
144 symbols.append(symbol.capitalize())
145 pos.append(list(map(float, [x, y, z])))
146 atoms = Atoms(symbols, np.array(pos) * Bohr)
147 continue
149 # Energy
150 ematch = _energy_re.match(line)
151 if ematch is not None:
152 energy = float(ematch.group(1)) * Hartree
154 # MPn energy. Supplants energy parsed above.
155 elif line.strip().startswith('TOTAL ENERGY'):
156 energy = float(line.strip().split()[-1]) * Hartree
158 # Higher-level energy (e.g. coupled cluster)
159 # Supplants energies parsed above.
160 elif line.strip().startswith('THE FOLLOWING METHOD AND ENERGY'):
161 energy = float(fd.readline().strip().split()[-1]) * Hartree
163 # Gradients
164 elif _grad_re.match(line):
165 for _ in range(3):
166 fd.readline()
167 grad = []
168 while True:
169 atom = _atom_re.match(fd.readline())
170 if atom is None:
171 break
172 grad.append(list(map(float, atom.groups()[2:])))
173 forces = -np.array(grad) * Hartree / Bohr
174 elif _dipole_re.match(line):
175 dipole = np.array(list(map(float, fd.readline().split()[:3])))
176 dipole *= Debye
178 atoms.calc = SinglePointCalculator(atoms, energy=energy,
179 forces=forces, dipole=dipole)
180 return atoms
183def read_gamess_us_punch(fd):
184 atoms = None
185 energy = None
186 forces = None
187 dipole = None
188 for line in fd:
189 if line.strip() == '$DATA':
190 symbols = []
191 pos = []
192 while line.strip() != '$END':
193 line = fd.readline()
194 atom = _atom_re.match(line)
195 if atom is None:
196 # The basis set specification is interlaced with the
197 # molecular geometry. We don't care about the basis
198 # set, so ignore lines that don't match the pattern.
199 continue
200 symbols.append(atom.group(1).capitalize())
201 pos.append(list(map(float, atom.group(3, 4, 5))))
202 atoms = Atoms(symbols, np.array(pos))
203 elif line.startswith('E('):
204 energy = float(line.split()[1][:-1]) * Hartree
205 elif line.strip().startswith('DIPOLE'):
206 dipole = np.array(list(map(float, line.split()[1:]))) * Debye
207 elif line.strip() == '$GRAD':
208 # The gradient block also contains the energy, which we prefer
209 # over the energy obtained above because it is more likely to
210 # be consistent with the gradients. It probably doesn't actually
211 # make a difference though.
212 energy = float(fd.readline().split()[1]) * Hartree
213 grad = []
214 while line.strip() != '$END':
215 line = fd.readline()
216 atom = _atom_re.match(line)
217 if atom is None:
218 continue
219 grad.append(list(map(float, atom.group(3, 4, 5))))
220 forces = -np.array(grad) * Hartree / Bohr
222 atoms.calc = SinglePointCalculator(atoms, energy=energy, forces=forces,
223 dipole=dipole)
225 return atoms
228def clean_userscr(userscr, prefix):
229 for fname in os.listdir(userscr):
230 tokens = fname.split('.')
231 if tokens[0] == prefix and tokens[-1] != 'bak':
232 fold = os.path.join(userscr, fname)
233 os.rename(fold, fold + '.bak')
236def get_userscr(prefix, command):
237 prefix_test = prefix + '_test'
238 command = command.replace('PREFIX', prefix_test)
239 with workdir(prefix_test, mkdir=True):
240 try:
241 call(command, shell=True, timeout=2)
242 except TimeoutExpired:
243 pass
245 try:
246 with open(prefix_test + '.log') as fd:
247 for line in fd:
248 if line.startswith('GAMESS supplementary output files'):
249 return ' '.join(line.split(' ')[8:]).strip()
250 except FileNotFoundError:
251 return None
253 return None