Coverage for /builds/kinetik161/ase/ase/io/gamess_us.py: 34.66%
176 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1import os
2import re
3from copy import deepcopy
4from subprocess import TimeoutExpired, call
6import numpy as np
8from ase import Atoms
9from ase.calculators.singlepoint import SinglePointCalculator
10from ase.units import Bohr, Debye, Hartree
11from ase.utils import workdir
14def _format_value(val):
15 if isinstance(val, bool):
16 return '.t.' if val else '.f.'
17 return str(val).upper()
20def _write_block(name, args):
21 out = [f' ${name.upper()}']
22 for key, val in args.items():
23 out.append(f' {key.upper()}={_format_value(val)}')
24 out.append(' $END')
25 return '\n'.join(out)
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)
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)
59_xc = dict(LDA='SVWN')
62def write_gamess_us_in(fd, atoms, properties=None, **params):
63 params = deepcopy(params)
65 if properties is None:
66 properties = ['energy']
68 # set RUNTYP from properties iff value not provided by the user
69 contrl = params.pop('contrl', {})
70 if 'runtyp' not in contrl:
71 if 'forces' in properties:
72 contrl['runtyp'] = 'gradient'
73 else:
74 contrl['runtyp'] = 'energy'
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())
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
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'
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'
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')
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))
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\)')
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
147 # Energy
148 ematch = _energy_re.match(line)
149 if ematch is not None:
150 energy = float(ematch.group(1)) * Hartree
152 # MPn energy. Supplants energy parsed above.
153 elif line.strip().startswith('TOTAL ENERGY'):
154 energy = float(line.strip().split()[-1]) * Hartree
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
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
176 atoms.calc = SinglePointCalculator(atoms, energy=energy,
177 forces=forces, dipole=dipole)
178 return atoms
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
220 atoms.calc = SinglePointCalculator(atoms, energy=energy, forces=forces,
221 dipole=dipole)
223 return atoms
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')
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
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
251 return None