Coverage for /builds/debichem-team/python-ase/ase/io/gaussian.py : 84.49%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import re
2import warnings
3from collections.abc import Iterable
4from copy import deepcopy
6import numpy as np
7from ase import Atoms
8from ase.calculators.calculator import InputError, Calculator
9from ase.calculators.gaussian import Gaussian
10from ase.calculators.singlepoint import SinglePointCalculator
11from ase.data import atomic_masses_iupac2016, chemical_symbols
12from ase.io import ParseError
13from ase.io.zmatrix import parse_zmatrix
14from ase.units import Bohr, Hartree
16_link0_keys = [
17 'mem',
18 'chk',
19 'oldchk',
20 'schk',
21 'rwf',
22 'oldmatrix',
23 'oldrawmatrix',
24 'int',
25 'd2e',
26 'save',
27 'nosave',
28 'errorsave',
29 'cpu',
30 'nprocshared',
31 'gpucpu',
32 'lindaworkers',
33 'usessh',
34 'ssh',
35 'debuglinda',
36]
39_link0_special = [
40 'kjob',
41 'subst',
42]
45# Certain problematic methods do not provide well-defined potential energy
46# surfaces, because these "composite" methods involve geometry optimization
47# and/or vibrational frequency analysis. In addition, the "energy" calculated
48# by these methods are typically ZPVE corrected and/or temperature dependent
49# free energies.
50_problem_methods = [
51 'cbs-4m', 'cbs-qb3', 'cbs-apno',
52 'g1', 'g2', 'g3', 'g4', 'g2mp2', 'g3mp2', 'g3b3', 'g3mp2b3', 'g4mp4',
53 'w1', 'w1u', 'w1bd', 'w1ro',
54]
57_xc_to_method = dict(
58 pbe='pbepbe',
59 pbe0='pbe1pbe',
60 hse06='hseh1pbe',
61 hse03='ohse2pbe',
62 lda='svwn', # gaussian "knows about" LSDA, but maybe not LDA.
63 tpss='tpsstpss',
64 revtpss='revtpssrevtpss',
65)
67_nuclear_prop_names = ['spin', 'zeff', 'qmom', 'nmagm', 'znuc',
68 'radnuclear', 'iso']
71def _get_molecule_spec(atoms, nuclear_props):
72 ''' Generate the molecule specification section to write
73 to the Gaussian input file, from the Atoms object and dict
74 of nuclear properties'''
75 molecule_spec = []
76 for i, atom in enumerate(atoms):
77 symbol_section = atom.symbol + '('
78 # Check whether any nuclear properties of the atom have been set,
79 # and if so, add them to the symbol section.
80 nuclear_props_set = False
81 for keyword, array in nuclear_props.items():
82 if array is not None and array[i] is not None:
83 string = keyword + '=' + str(array[i]) + ', '
84 symbol_section += string
85 nuclear_props_set = True
87 # Check whether the mass of the atom has been modified,
88 # and if so, add it to the symbol section:
89 mass_set = False
90 symbol = atom.symbol
91 expected_mass = atomic_masses_iupac2016[chemical_symbols.index(
92 symbol)]
93 if expected_mass != atoms[i].mass:
94 mass_set = True
95 string = 'iso' + '=' + str(atoms[i].mass)
96 symbol_section += string
98 if nuclear_props_set or mass_set:
99 symbol_section = symbol_section.strip(', ')
100 symbol_section += ')'
101 else:
102 symbol_section = symbol_section.strip('(')
104 # Then attach the properties appropriately
105 # this formatting was chosen for backwards compatibility reasons, but
106 # it would probably be better to
107 # 1) Ensure proper spacing between entries with explicit spaces
108 # 2) Use fewer columns for the element
109 # 3) Use 'e' (scientific notation) instead of 'f' for positions
110 molecule_spec.append('{:<10s}{:20.10f}{:20.10f}{:20.10f}'.format(
111 symbol_section, *atom.position))
113 # unit cell vectors, in case of periodic boundary conditions
114 for ipbc, tv in zip(atoms.pbc, atoms.cell):
115 if ipbc:
116 molecule_spec.append('TV {:20.10f}{:20.10f}{:20.10f}'.format(*tv))
118 molecule_spec.append('')
120 return molecule_spec
123def _format_output_type(output_type):
124 ''' Given a letter: output_type, return
125 a string formatted for a gaussian input file'''
126 # Change output type to P if it has been set to T, because ASE
127 # does not support reading info from 'terse' output files.
128 if output_type is None or output_type == '' or 't' in output_type.lower():
129 output_type = 'P'
131 return '#{}'.format(output_type)
134def _check_problem_methods(method):
135 ''' Check method string for problem methods and warn appropriately'''
136 if method.lower() in _problem_methods:
137 warnings.warn(
138 'The requested method, {}, is a composite method. Composite '
139 'methods do not have well-defined potential energy surfaces, '
140 'so the energies, forces, and other properties returned by '
141 'ASE may not be meaningful, or they may correspond to a '
142 'different geometry than the one provided. '
143 'Please use these methods with caution.'.format(method)
144 )
147def _pop_link0_params(params):
148 '''Takes the params dict and returns a dict with the link0 keywords
149 removed, and a list containing the link0 keywords and options
150 to be written to the gaussian input file'''
151 params = deepcopy(params)
152 out = []
153 # Remove keywords from params, so they are not set again later in
154 # route section
155 for key in _link0_keys:
156 if key not in params:
157 continue
158 val = params.pop(key)
159 if not val or (isinstance(val, str) and key.lower() == val.lower()):
160 out.append('%{}'.format(key))
161 else:
162 out.append('%{}={}'.format(key, val))
164 # These link0 keywords have a slightly different syntax
165 for key in _link0_special:
166 if key not in params:
167 continue
168 val = params.pop(key)
169 if not isinstance(val, str) and isinstance(val, Iterable):
170 val = ' '.join(val)
171 out.append('%{} L{}'.format(key, val))
173 return params, out
176def _format_method_basis(output_type, method, basis, fitting_basis):
177 output_string = ""
178 if basis and method and fitting_basis:
179 output_string = '{} {}/{}/{} ! ASE formatted method and basis'.format(
180 output_type, method, basis, fitting_basis)
181 elif basis and method:
182 output_string = '{} {}/{} ! ASE formatted method and basis'.format(
183 output_type, method, basis)
184 else:
185 output_string = '{}'.format(output_type)
186 for value in [method, basis]:
187 if value is not None:
188 output_string += ' {}'.format(value)
189 return output_string
192def _format_route_params(params):
193 '''Get keywords and values from the params dictionary and return
194 as a list of lines to add to the gaussian input file'''
195 out = []
196 for key, val in params.items():
197 # assume bare keyword if val is falsey, i.e. '', None, False, etc.
198 # also, for backwards compatibility: assume bare keyword if key and
199 # val are the same
200 if not val or (isinstance(val, str) and key.lower() == val.lower()):
201 out.append(key)
202 elif not isinstance(val, str) and isinstance(val, Iterable):
203 out.append('{}({})'.format(key, ','.join(val)))
204 else:
205 out.append('{}({})'.format(key, val))
206 return out
209def _format_addsec(addsec):
210 '''Format addsec string as a list of lines to be added to the gaussian
211 input file.'''
212 out = []
213 if addsec is not None:
214 out.append('')
215 if isinstance(addsec, str):
216 out.append(addsec)
217 elif isinstance(addsec, Iterable):
218 out += list(addsec)
219 return out
222def _format_basis_set(basis, basisfile, basis_set):
223 '''Format either: the basis set filename (basisfile), the basis set file
224 contents (from reading basisfile), or the basis_set text as a list of
225 strings to be added to the gaussian input file.'''
226 out = []
227 # if basis='gen', set basisfile. Either give a path to a basisfile, or
228 # read in the provided file and paste it verbatim
229 if basisfile is not None:
230 if basisfile[0] == '@':
231 out.append(basisfile)
232 else:
233 with open(basisfile, 'r') as fd:
234 out.append(fd.read())
235 elif basis_set is not None:
236 out.append(basis_set)
237 else:
238 if basis is not None and basis.lower() == 'gen':
239 raise InputError('Please set basisfile or basis_set')
240 return out
243def write_gaussian_in(fd, atoms, properties=['energy'],
244 method=None, basis=None, fitting_basis=None,
245 output_type='P', basisfile=None, basis_set=None,
246 xc=None, charge=None, mult=None, extra=None,
247 ioplist=None, addsec=None, spinlist=None,
248 zefflist=None, qmomlist=None, nmagmlist=None,
249 znuclist=None, radnuclearlist=None,
250 **params):
251 '''
252 Generates a Gaussian input file
254 Parameters
255 -----------
256 fd: file-like
257 where the Gaussian input file will be written
258 atoms: Atoms
259 Structure to write to the input file
260 properties: list
261 Properties to calculate
262 method: str
263 Level of theory to use, e.g. ``hf``, ``ccsd``, ``mp2``, or ``b3lyp``.
264 Overrides ``xc`` (see below).
265 xc: str
266 Level of theory to use. Translates several XC functionals from
267 their common name (e.g. ``PBE``) to their internal Gaussian name
268 (e.g. ``PBEPBE``).
269 basis: str
270 The basis set to use. If not provided, no basis set will be requested,
271 which usually results in ``STO-3G``. Maybe omitted if basisfile is set
272 (see below).
273 fitting_basis: str
274 The name of the fitting basis set to use.
275 output_type: str
276 Level of output to record in the Gaussian
277 output file - this may be ``N``- normal or ``P`` -
278 additional.
279 basisfile: str
280 The name of the basis file to use. If a value is provided, basis may
281 be omitted (it will be automatically set to 'gen')
282 basis_set: str
283 The basis set definition to use. This is an alternative
284 to basisfile, and would be the same as the contents
285 of such a file.
286 charge: int
287 The system charge. If not provided, it will be automatically
288 determined from the ``Atoms`` object’s initial_charges.
289 mult: int
290 The system multiplicity (``spin + 1``). If not provided, it will be
291 automatically determined from the ``Atoms`` object’s
292 ``initial_magnetic_moments``.
293 extra: str
294 Extra lines to be included in the route section verbatim.
295 It should not be necessary to use this, but it is included for
296 backwards compatibility.
297 ioplist: list
298 A collection of IOPs definitions to be included in the route line.
299 addsec: str
300 Text to be added after the molecular geometry specification, e.g. for
301 defining masses with ``freq=ReadIso``.
302 spinlist: list
303 A list of nuclear spins to be added into the nuclear
304 propeties section of the molecule specification.
305 zefflist: list
306 A list of effective charges to be added into the nuclear
307 propeties section of the molecule specification.
308 qmomlist: list
309 A list of nuclear quadropole moments to be added into
310 the nuclear propeties section of the molecule
311 specification.
312 nmagmlist: list
313 A list of nuclear magnetic moments to be added into
314 the nuclear propeties section of the molecule
315 specification.
316 znuclist: list
317 A list of nuclear charges to be added into the nuclear
318 propeties section of the molecule specification.
319 radnuclearlist: list
320 A list of nuclear radii to be added into the nuclear
321 propeties section of the molecule specification.
322 params: dict
323 Contains any extra keywords and values that will be included in either
324 the link0 section or route section of the gaussian input file.
325 To be included in the link0 section, the keyword must be one of the
326 following: ``mem``, ``chk``, ``oldchk``, ``schk``, ``rwf``,
327 ``oldmatrix``, ``oldrawmatrix``, ``int``, ``d2e``, ``save``,
328 ``nosave``, ``errorsave``, ``cpu``, ``nprocshared``, ``gpucpu``,
329 ``lindaworkers``, ``usessh``, ``ssh``, ``debuglinda``.
330 Any other keywords will be placed (along with their values) in the
331 route section.
332 '''
334 params = deepcopy(params)
336 if properties is None:
337 properties = ['energy']
339 output_type = _format_output_type(output_type)
341 # basis can be omitted if basisfile is provided
342 if basis is None:
343 if basisfile is not None or basis_set is not None:
344 basis = 'gen'
346 # determine method from xc if it is provided
347 if method is None:
348 if xc is not None:
349 method = _xc_to_method.get(xc.lower(), xc)
351 # If the user requests a problematic method, rather than raising an error
352 # or proceeding blindly, give the user a warning that the results parsed
353 # by ASE may not be meaningful.
354 if method is not None:
355 _check_problem_methods(method)
357 # determine charge from initial charges if not passed explicitly
358 if charge is None:
359 charge = atoms.get_initial_charges().sum()
361 # determine multiplicity from initial magnetic moments
362 # if not passed explicitly
363 if mult is None:
364 mult = atoms.get_initial_magnetic_moments().sum() + 1
366 # set up link0 arguments
367 out = []
368 params, link0_list = _pop_link0_params(params)
369 out.extend(link0_list)
371 # begin route line
372 # note: unlike in old calculator, each route keyword is put on its own
373 # line.
374 out.append(_format_method_basis(output_type, method, basis, fitting_basis))
376 # If the calculator's parameter dictionary contains an isolist, we ignore
377 # this - it is up to the user to attach this info as the atoms' masses
378 # if they wish for it to be used:
379 params.pop('isolist', None)
381 # Any params left will belong in the route section of the file:
382 out.extend(_format_route_params(params))
384 if ioplist is not None:
385 out.append('IOP(' + ', '.join(ioplist) + ')')
387 # raw list of explicit keywords for backwards compatibility
388 if extra is not None:
389 out.append(extra)
391 # Add 'force' iff the user requested forces, since Gaussian crashes when
392 # 'force' is combined with certain other keywords such as opt and irc.
393 if 'forces' in properties and 'force' not in params:
394 out.append('force')
396 # header, charge, and mult
397 out += ['', 'Gaussian input prepared by ASE', '',
398 '{:.0f} {:.0f}'.format(charge, mult)]
400 # make dict of nuclear properties:
401 nuclear_props = {'spin': spinlist, 'zeff': zefflist, 'qmom': qmomlist,
402 'nmagm': nmagmlist, 'znuc': znuclist,
403 'radnuclear': radnuclearlist}
404 nuclear_props = {k: v for k, v in nuclear_props.items() if v is not None}
406 # atomic positions and nuclear properties:
407 molecule_spec = _get_molecule_spec(atoms, nuclear_props)
408 for line in molecule_spec:
409 out.append(line)
411 out.extend(_format_basis_set(basis, basisfile, basis_set))
413 out.extend(_format_addsec(addsec))
415 out += ['', '']
416 fd.write('\n'.join(out))
419# Regexp for reading an input file:
421_re_link0 = re.compile(r'^\s*%([^\=\)\(!]+)=?([^\=\)\(!]+)?(!.+)?')
422# Link0 lines are in the format:
423# '% keyword = value' or '% keyword'
424# (with or without whitespaces)
426_re_output_type = re.compile(r'^\s*#\s*([NPTnpt]?)\s*')
427# The start of the route section begins with a '#', and then may
428# be followed by the desired level of output in the output file: P, N or T.
430_re_method_basis = re.compile(
431 r"\s*([\w-]+)\s*\/([^/=!]+)([\/]([^!]+))?\s*(!.+)?")
432# Matches method, basis and optional fitting basis in the format:
433# method/basis/fitting_basis ! comment
434# They will appear in this format if the Gaussian file has been generated
435# by ASE using a calculator with the basis and method keywords set.
437_re_chgmult = re.compile(r'^\s*[+-]?\d+(?:,\s*|\s+)[+-]?\d+\s*$')
438# This is a bit more complex of a regex than we typically want, but it
439# can be difficult to determine whether a line contains the charge and
440# multiplicity, rather than just another route keyword. By making sure
441# that the line contains exactly two *integers*, separated by either
442# a comma (and possibly whitespace) or some amount of whitespace, we
443# can be more confident that we've actually found the charge and multiplicity.
445_re_nuclear_props = re.compile(r'\(([^\)]+)\)')
446# Matches the nuclear properties, which are contained in parantheses.
448# The following methods are used in GaussianConfiguration's
449# parse_gaussian_input method:
452def _get_link0_param(link0_match):
453 '''Gets link0 keyword and option from a re.Match and returns them
454 in a dictionary format'''
455 value = link0_match.group(2)
456 if value is not None:
457 value = value.strip()
458 else:
459 value = ''
460 return {link0_match.group(1).lower().strip(): value.lower()}
463def _get_all_link0_params(link0_section):
464 ''' Given a string link0_section which contains the link0
465 section of a gaussian input file, returns a dictionary of
466 keywords and values'''
467 parameters = {}
468 for line in link0_section:
469 link0_match = _re_link0.match(line)
470 link0_param = _get_link0_param(link0_match)
471 if link0_param is not None:
472 parameters.update(link0_param)
473 return parameters
476def _convert_to_symbol(string):
477 '''Converts an input string into a format
478 that can be input to the 'symbol' parameter of an
479 ASE Atom object (can be a chemical symbol (str)
480 or an atomic number (int)).
481 This is achieved by either stripping any
482 integers from the string, or converting a string
483 containing an atomic number to integer type'''
484 symbol = _validate_symbol_string(string)
485 if symbol.isnumeric():
486 atomic_number = int(symbol)
487 symbol = chemical_symbols[atomic_number]
488 else:
489 match = re.match(r'([A-Za-z]+)', symbol)
490 symbol = match.group(1)
491 return symbol
494def _validate_symbol_string(string):
495 if "-" in string:
496 raise ParseError("ERROR: Could not read the Gaussian input file, as"
497 " molecule specifications for molecular mechanics "
498 "calculations are not supported.")
499 return string
502def _get_key_value_pairs(line):
503 '''Reads a line of a gaussian input file, which contains keywords and options
504 separated according to the rules of the route section.
506 Parameters
507 ----------
508 line (string)
509 A line of a gaussian input file.
511 Returns
512 ---------
513 params (dict)
514 Contains the keywords and options found in the line.
515 '''
516 params = {}
517 line = line.strip(' #')
518 line = line.split('!')[0] # removes any comments
519 # First, get the keywords and options sepatated with
520 # parantheses:
521 match_iterator = re.finditer(r'\(([^\)]+)\)', line)
522 index_ranges = []
523 for match in match_iterator:
524 index_range = [match.start(0), match.end(0)]
525 options = match.group(1)
526 # keyword is last word in previous substring:
527 keyword_string = line[:match.start(0)]
528 keyword_match_iter = [k for k in re.finditer(
529 r'[^\,/\s]+', keyword_string) if k.group() != '=']
530 keyword = keyword_match_iter[-1].group().strip(' =')
531 index_range[0] = keyword_match_iter[-1].start()
532 params.update({keyword.lower(): options.lower()})
533 index_ranges.append(index_range)
535 # remove from the line the keywords and options that we have saved:
536 index_ranges.reverse()
537 for index_range in index_ranges:
538 start = index_range[0]
539 stop = index_range[1]
540 line = line[0: start:] + line[stop + 1::]
542 # Next, get the keywords and options separated with
543 # an equals sign, and those without an equals sign
544 # must be keywords without options:
546 # remove any whitespaces around '=':
547 line = re.sub(r'\s*=\s*', '=', line)
548 line = [x for x in re.split(r'[\s,\/]', line) if x != '']
550 for s in line:
551 if '=' in s:
552 s = s.split('=')
553 keyword = s.pop(0)
554 options = s.pop(0)
555 params.update({keyword.lower(): options.lower()})
556 else:
557 if len(s) > 0:
558 params.update({s.lower(): None})
560 return params
563def _get_route_params(line):
564 '''Reads keywords and values from a line in
565 a Gaussian input file's route section,
566 and returns them as a dictionary'''
567 method_basis_match = _re_method_basis.match(line)
568 if method_basis_match:
569 params = {}
570 ase_gen_comment = '! ASE formatted method and basis'
571 if method_basis_match.group(5) == ase_gen_comment:
572 params['method'] = method_basis_match.group(1).strip().lower()
573 params['basis'] = method_basis_match.group(2).strip().lower()
574 if method_basis_match.group(4):
575 params['fitting_basis'] = method_basis_match.group(
576 4).strip().lower()
577 return params
579 return _get_key_value_pairs(line)
582def _get_all_route_params(route_section):
583 ''' Given a string: route_section which contains the route
584 section of a gaussian input file, returns a dictionary of
585 keywords and values'''
586 parameters = {}
588 for line in route_section:
589 output_type_match = _re_output_type.match(line)
590 if not parameters.get('output_type') and output_type_match:
591 line = line.strip(output_type_match.group(0))
592 parameters.update(
593 {'output_type': output_type_match.group(1).lower()})
594 # read route section
595 route_params = _get_route_params(line)
596 if route_params is not None:
597 parameters.update(route_params)
598 return parameters
601def _get_charge_mult(chgmult_section):
602 '''return a dict with the charge and multiplicity from
603 a list chgmult_section that contains the charge and multiplicity
604 line, read from a gaussian input file'''
605 chgmult_match = _re_chgmult.match(str(chgmult_section))
606 try:
607 chgmult = chgmult_match.group(0).split()
608 return {'charge': int(chgmult[0]), 'mult': int(chgmult[1])}
609 except (IndexError, AttributeError):
610 raise ParseError("ERROR: Could not read the charge and multiplicity "
611 "from the Gaussian input file. These must be 2 "
612 "integers separated with whitespace or a comma.")
615def _get_nuclear_props(line):
616 ''' Reads any info in parantheses in the line and returns
617 a dictionary of the nuclear properties.'''
618 nuclear_props_match = _re_nuclear_props.search(line)
619 nuclear_props = {}
620 if nuclear_props_match:
621 nuclear_props = _get_key_value_pairs(nuclear_props_match.group(1))
622 updated_nuclear_props = {}
623 for key, value in nuclear_props.items():
624 if value.isnumeric():
625 value = int(value)
626 else:
627 value = float(value)
628 if key not in _nuclear_prop_names:
629 if "fragment" in key:
630 warnings.warn("Fragments are not "
631 "currently supported.")
632 warnings.warn("The following nuclear properties "
633 "could not be saved: {}".format(
634 {key: value}))
635 else:
636 updated_nuclear_props[key] = value
637 nuclear_props = updated_nuclear_props
639 for k in _nuclear_prop_names:
640 if k not in nuclear_props:
641 nuclear_props[k] = None
643 return nuclear_props
646def _get_atoms_info(line):
647 '''Returns the symbol and position of an atom from a line
648 in the molecule specification section'''
649 nuclear_props_match = _re_nuclear_props.search(line)
650 if nuclear_props_match:
651 line = line.replace(nuclear_props_match.group(0), '')
652 tokens = line.split()
653 symbol = _convert_to_symbol(tokens[0])
654 pos = list(tokens[1:])
656 return symbol, pos
659def _get_cartesian_atom_coords(symbol, pos):
660 '''Returns the coordinates: pos as a list of floats if they
661 are cartesian, and not in z-matrix format'''
662 if len(pos) < 3 or (pos[0] == '0' and symbol != 'TV'):
663 # In this case, we have a z-matrix definition, so
664 # no cartesian coords.
665 return
666 elif len(pos) > 3:
667 raise ParseError("ERROR: Gaussian input file could "
668 "not be read as freeze codes are not"
669 " supported. If using cartesian "
670 "coordinates, these must be "
671 "given as 3 numbers separated "
672 "by whitespace.")
673 else:
674 try:
675 return list(map(float, pos))
676 except ValueError:
677 raise(ParseError(
678 "ERROR: Molecule specification in"
679 "Gaussian input file could not be read"))
682def _get_zmatrix_line(line):
683 ''' Converts line into the format needed for it to
684 be added to the z-matrix contents '''
685 line_list = line.split()
686 if len(line_list) == 8 and line_list[7] == '1':
687 raise ParseError(
688 "ERROR: Could not read the Gaussian input file"
689 ", as the alternative Z-matrix format using "
690 "two bond angles instead of a bond angle and "
691 "a dihedral angle is not supported.")
692 return(line.strip() + '\n')
695def _read_zmatrix(zmatrix_contents, zmatrix_vars=None):
696 ''' Reads a z-matrix (zmatrix_contents) using its list of variables
697 (zmatrix_vars), and returns atom positions and symbols '''
698 try:
699 atoms = parse_zmatrix(zmatrix_contents, defs=zmatrix_vars)
700 except (ValueError, AssertionError) as e:
701 raise ParseError("Failed to read Z-matrix from "
702 "Gaussian input file: ", e)
703 except KeyError as e:
704 raise ParseError("Failed to read Z-matrix from "
705 "Gaussian input file, as symbol: {}"
706 "could not be recognised. Please make "
707 "sure you use element symbols, not "
708 "atomic numbers in the element labels.".format(e))
709 positions = atoms.positions
710 symbols = atoms.get_chemical_symbols()
711 return positions, symbols
714def _get_nuclear_props_for_all_atoms(nuclear_props):
715 ''' Returns the nuclear properties for all atoms as a dictionary,
716 in the format needed for it to be added to the parameters dictionary.'''
717 params = {k + 'list': [] for k in _nuclear_prop_names}
718 for dictionary in nuclear_props:
719 for key, value in dictionary.items():
720 params[key + 'list'].append(value)
722 for key, array in params.items():
723 values_set = False
724 for value in array:
725 if value is not None:
726 values_set = True
727 if not values_set:
728 params[key] = None
729 return params
732def _get_atoms_from_molspec(molspec_section):
733 ''' Takes a string: molspec_section which contains the molecule
734 specification section of a gaussian input file, and returns an atoms
735 object that represents this.'''
736 # These will contain info that will be attached to the Atoms object:
737 symbols = []
738 positions = []
739 pbc = np.zeros(3, dtype=bool)
740 cell = np.zeros((3, 3))
741 npbc = 0
743 # Will contain a dictionary of nuclear properties for each atom,
744 # that will later be saved to the parameters dict:
745 nuclear_props = []
747 # Info relating to the z-matrix definition (if set)
748 zmatrix_type = False
749 zmatrix_contents = ""
750 zmatrix_var_section = False
751 zmatrix_vars = ""
753 for line in molspec_section:
754 # Remove any comments and replace '/' and ',' with whitespace,
755 # as these are equivalent:
756 line = line.split('!')[0].replace('/', ' ').replace(',', ' ')
757 if (line.split()):
758 if zmatrix_type:
759 # Save any variables set when defining the z-matrix:
760 if zmatrix_var_section:
761 zmatrix_vars += line.strip() + '\n'
762 continue
763 elif 'variables' in line.lower():
764 zmatrix_var_section = True
765 continue
766 elif 'constants' in line.lower():
767 zmatrix_var_section = True
768 warnings.warn("Constants in the optimisation are "
769 "not currently supported. Instead "
770 "setting constants as variables.")
771 continue
773 symbol, pos = _get_atoms_info(line)
774 current_nuclear_props = _get_nuclear_props(line)
776 if not zmatrix_type:
777 pos = _get_cartesian_atom_coords(symbol, pos)
778 if pos is None:
779 zmatrix_type = True
781 if symbol.upper() == 'TV' and pos is not None:
782 pbc[npbc] = True
783 cell[npbc] = pos
784 npbc += 1
785 else:
786 nuclear_props.append(current_nuclear_props)
787 if not zmatrix_type:
788 symbols.append(symbol)
789 positions.append(pos)
791 if zmatrix_type:
792 zmatrix_contents += _get_zmatrix_line(line)
794 # Now that we are past the molecule spec. section, we can read
795 # the entire z-matrix (if set):
796 if len(positions) == 0:
797 if zmatrix_type:
798 if zmatrix_vars == '':
799 zmatrix_vars = None
800 positions, symbols = _read_zmatrix(
801 zmatrix_contents, zmatrix_vars)
803 try:
804 atoms = Atoms(symbols, positions, pbc=pbc, cell=cell)
805 except (IndexError, ValueError, KeyError) as e:
806 raise ParseError("ERROR: Could not read the Gaussian input file, "
807 "due to a problem with the molecule "
808 "specification: {}".format(e))
810 nuclear_props = _get_nuclear_props_for_all_atoms(nuclear_props)
812 return atoms, nuclear_props
815def _get_readiso_param(parameters):
816 ''' Returns a tuple containing the frequency
817 keyword and its options, if the frequency keyword is
818 present in parameters and ReadIso is one of its options'''
819 freq_options = parameters.get('freq', None)
820 if freq_options:
821 freq_name = 'freq'
822 else:
823 freq_options = parameters.get('frequency', None)
824 freq_name = 'frequency'
825 if freq_options is not None:
826 if ('readiso' or 'readisotopes') in freq_options:
827 return freq_name, freq_options
828 return None, None
831def _get_readiso_info(line, parameters):
832 '''Reads the temperature, pressure and scale from the first line
833 of a ReadIso section of a Gaussian input file. Returns these in
834 a dictionary.'''
835 readiso_params = {}
836 if _get_readiso_param(parameters)[0] is not None:
837 # when count_iso is 0 we are in the line where
838 # temperature, pressure, [scale] is saved
839 line = line.replace(
840 '[', '').replace(']', '')
841 tokens = line.strip().split()
842 try:
843 readiso_params['temperature'] = tokens[0]
844 readiso_params['pressure'] = tokens[1]
845 readiso_params['scale'] = tokens[2]
846 except IndexError:
847 pass
848 if readiso_params != {}:
850 return readiso_params
853def _delete_readiso_param(parameters):
854 '''Removes the readiso parameter from the parameters dict'''
855 parameters = deepcopy(parameters)
856 freq_name, freq_options = _get_readiso_param(parameters)
857 if freq_name is not None:
858 if 'readisotopes' in freq_options:
859 iso_name = 'readisotopes'
860 else:
861 iso_name = 'readiso'
862 freq_options = [v.group() for v in re.finditer(
863 r'[^\,/\s]+', freq_options)]
864 freq_options.remove(iso_name)
865 new_freq_options = ''
866 for v in freq_options:
867 new_freq_options += v + ' '
868 if new_freq_options == '':
869 new_freq_options = None
870 else:
871 new_freq_options = new_freq_options.strip()
872 parameters[freq_name] = new_freq_options
873 return parameters
876def _update_readiso_params(parameters, symbols):
877 ''' Deletes the ReadIso option from the route section as we
878 write out the masses in the nuclear properties section
879 instead of the ReadIso section.
880 Ensures the masses array is the same length as the
881 symbols array. This is necessary due to the way the
882 ReadIso section is defined:
883 The mass of each atom is listed on a separate line, in
884 order of appearance in the molecule spec. A blank line
885 indicates not to modify the mass for that atom.
886 But you do not have to leave blank lines equal to the
887 remaining atoms after you finsihed setting masses.
888 E.g. if you had 10 masses and only want to set the mass
889 for the first atom, you don't have to leave 9 blank lines
890 after it.
891 '''
892 parameters = _delete_readiso_param(parameters)
893 if parameters.get('isolist') is not None:
894 if len(parameters['isolist']) < len(symbols):
895 for i in range(0, len(symbols) - len(parameters['isolist'])):
896 parameters['isolist'].append(None)
897 if all(m is None for m in parameters['isolist']):
898 parameters['isolist'] = None
900 return parameters
903def _validate_params(parameters):
904 '''Checks whether all of the required parameters exist in the
905 parameters dict and whether it contains any unsupported settings
906 '''
907 # Check for unsupported settings
908 unsupported_settings = {
909 "z-matrix", "modredun", "modredundant", "addredundant", "addredun",
910 "readopt", "rdopt"}
912 for s in unsupported_settings:
913 for v in parameters.values():
914 if v is not None and s in str(v):
915 raise ParseError(
916 "ERROR: Could not read the Gaussian input file"
917 ", as the option: {} is currently unsupported."
918 .format(s))
920 for k in list(parameters.keys()):
921 if "popt" in k:
922 parameters["opt"] = parameters.pop(k)
923 warnings.warn("The option {} is currently unsupported. "
924 "This has been replaced with {}."
925 .format("POpt", "opt"))
926 return
929def _get_extra_section_params(extra_section, parameters, atoms):
930 ''' Takes a list of strings: extra_section, which contains
931 the 'extra' lines in a gaussian input file. Also takes the parameters
932 that have been read so far, and the atoms that have been read from the
933 file.
934 Returns an updated version of the parameters dict, containing details from
935 this extra section. This may include the basis set definition or filename,
936 and/or the readiso section.'''
938 new_parameters = deepcopy(parameters)
939 # Will store the basis set definition (if set):
940 basis_set = ""
941 # This will indicate whether we have a readiso section:
942 readiso = _get_readiso_param(new_parameters)[0]
943 # This will indicate which line of the readiso section we are reading:
944 count_iso = 0
945 readiso_masses = []
947 for line in extra_section:
948 if line.split():
949 # check that the line isn't just a comment
950 if line.split()[0] == '!':
951 continue
952 line = line.strip().split('!')[0]
954 if len(line) > 0 and line[0] == '@':
955 # If the name of a basis file is specified, this line
956 # begins with a '@'
957 new_parameters['basisfile'] = line
958 elif readiso and count_iso < len(atoms.symbols) + 1:
959 # The ReadIso section has 1 more line than the number of
960 # symbols
961 if count_iso == 0 and line != '\n':
962 # The first line in the readiso section contains the
963 # temperature, pressure, scale. Here we save this:
964 readiso_info = _get_readiso_info(line, new_parameters)
965 if readiso_info is not None:
966 new_parameters.update(readiso_info)
967 # If the atom masses were set in the nuclear properties
968 # section, they will be overwritten by the ReadIso
969 # section
970 readiso_masses = []
971 count_iso += 1
972 elif count_iso > 0:
973 # The remaining lines in the ReadIso section are
974 # the masses of the atoms
975 try:
976 readiso_masses.append(float(line))
977 except ValueError:
978 readiso_masses.append(None)
979 count_iso += 1
980 else:
981 # If the rest of the file is not the ReadIso section,
982 # then it must be the definition of the basis set.
983 if new_parameters.get('basis', '') == 'gen' \
984 or 'gen' in new_parameters:
985 if line.strip() != '':
986 basis_set += line + '\n'
988 if readiso:
989 new_parameters['isolist'] = readiso_masses
990 new_parameters = _update_readiso_params(new_parameters, atoms.symbols)
992 # Saves the basis set definition to the parameters array if
993 # it has been set:
994 if basis_set != '':
995 new_parameters['basis_set'] = basis_set
996 new_parameters['basis'] = 'gen'
997 new_parameters.pop('gen', None)
999 return new_parameters
1002def _get_gaussian_in_sections(fd):
1003 ''' Reads a gaussian input file and returns
1004 a dictionary of the sections of the file - each dictionary
1005 value is a list of strings - each string is a line in that
1006 section. '''
1007 # These variables indicate which section of the
1008 # input file is currently being read:
1009 route_section = False
1010 atoms_section = False
1011 atoms_saved = False
1012 # Here we will store the sections of the file in a dictionary,
1013 # as lists of strings
1014 gaussian_sections = {'link0': [], 'route': [],
1015 'charge_mult': [], 'mol_spec': [], 'extra': []}
1017 for line in fd:
1018 line = line.strip(' ')
1019 link0_match = _re_link0.match(line)
1020 output_type_match = _re_output_type.match(line)
1021 chgmult_match = _re_chgmult.match(line)
1022 if link0_match:
1023 gaussian_sections['link0'].append(line)
1024 # The first blank line appears at the end of the route section
1025 # and a blank line appears at the end of the atoms section:
1026 elif line == '\n' and (route_section or atoms_section):
1027 route_section = False
1028 atoms_section = False
1029 elif output_type_match or route_section:
1030 route_section = True
1031 gaussian_sections['route'].append(line)
1032 elif chgmult_match:
1033 gaussian_sections['charge_mult'] = line
1034 # After the charge and multiplicty have been set, the
1035 # molecule specification section of the input file begins:
1036 atoms_section = True
1037 elif atoms_section:
1038 gaussian_sections['mol_spec'].append(line)
1039 atoms_saved = True
1040 elif atoms_saved:
1041 # Next we read the other sections below the molecule spec.
1042 # This may include the ReadIso section and the basis set
1043 # definition or filename
1044 gaussian_sections['extra'].append(line)
1046 return gaussian_sections
1049class GaussianConfiguration:
1051 def __init__(self, atoms, parameters):
1052 self.atoms = atoms.copy()
1053 self.parameters = deepcopy(parameters)
1055 def get_atoms(self):
1056 return self.atoms
1058 def get_parameters(self):
1059 return self.parameters
1061 def get_calculator(self):
1062 # Explicitly set parameters that must for security reasons not be
1063 # taken from file:
1064 calc = Gaussian(atoms=self.atoms, command=None, restart=None,
1065 ignore_bad_restart_file=Calculator._deprecated,
1066 label='Gaussian', directory='.', **self.parameters)
1067 return calc
1069 @ staticmethod
1070 def parse_gaussian_input(fd):
1071 '''Reads a gaussian input file into an atoms object and
1072 parameters dictionary.
1074 Parameters
1075 ----------
1076 fd: file-like
1077 Contains the contents of a gaussian input file
1079 Returns
1080 ---------
1081 GaussianConfiguration
1082 Contains an atoms object created using the structural
1083 information from the input file.
1084 Contains a parameters dictionary, which stores any
1085 keywords and options found in the link-0 and route
1086 sections of the input file.
1087 '''
1088 # The parameters dict to attach to the calculator
1089 parameters = {}
1091 file_sections = _get_gaussian_in_sections(fd)
1093 # Update the parameters dictionary with the keywords and values
1094 # from each section of the input file:
1095 parameters.update(_get_all_link0_params(file_sections['link0']))
1096 parameters.update(_get_all_route_params(file_sections['route']))
1097 parameters.update(_get_charge_mult(file_sections['charge_mult']))
1098 atoms, nuclear_props = _get_atoms_from_molspec(
1099 file_sections['mol_spec'])
1100 parameters.update(nuclear_props)
1101 parameters.update(_get_extra_section_params(
1102 file_sections['extra'], parameters, atoms))
1104 _validate_params(parameters)
1106 return GaussianConfiguration(atoms, parameters)
1109def read_gaussian_in(fd, attach_calculator=False):
1110 '''
1111 Reads a gaussian input file and returns an Atoms object.
1113 Parameters
1114 ----------
1115 fd: file-like
1116 Contains the contents of a gaussian input file
1118 attach_calculator: bool
1119 When set to ``True``, a Gaussian calculator will be
1120 attached to the Atoms object which is returned.
1121 This will mean that additional information is read
1122 from the input file (see below for details).
1124 Returns
1125 ----------
1126 Atoms
1127 An Atoms object containing the following information that has been
1128 read from the input file: symbols, positions, cell.
1130 Notes
1131 ----------
1132 Gaussian input files can be read where the atoms' locations are set with
1133 cartesian coordinates or as a z-matrix. Variables may be used in the
1134 z-matrix definition, but currently setting constants for constraining
1135 a geometry optimisation is not supported. Note that the `alternative`
1136 z-matrix definition where two bond angles are set instead of a bond angle
1137 and a dihedral angle is not currently supported.
1139 If the parameter ``attach_calculator`` is set to ``True``, then the Atoms
1140 object is returned with a Gaussian calculator attached.
1142 This Gaussian calculator will contain a parameters dictionary which will
1143 contain the Link 0 commands and the options and keywords set in the route
1144 section of the Gaussian input file, as well as:
1146 • The charge and multiplicity
1148 • The selected level of output
1150 • The method, basis set and (optional) fitting basis set.
1152 • Basis file name/definition if set
1154 • Nuclear properties
1156 • Masses set in the nuclear properties section or in the ReadIsotopes
1157 section (if ``freq=ReadIso`` is set). These will be saved in an array
1158 with keyword ``isolist``, in the parameters dictionary. For these to
1159 actually be used in calculations and/or written out to input files,
1160 the Atoms object's masses must be manually updated to these values
1161 (this is not done automatically)
1163 If the Gaussian input file has been generated by ASE's
1164 ``write_gaussian_in`` method, then the basis set, method and fitting
1165 basis will be saved under the ``basis``, ``method`` and ``fitting_basis``
1166 keywords in the parameters dictionary. Otherwise they are saved as
1167 keywords in the parameters dict.
1169 Currently this does not support reading of any other sections which would
1170 be found below the molecule specification that have not been mentioned
1171 here (such as the ModRedundant section).
1172 '''
1173 gaussian_input = GaussianConfiguration.parse_gaussian_input(fd)
1174 atoms = gaussian_input.get_atoms()
1176 if attach_calculator:
1177 atoms.calc = gaussian_input.get_calculator()
1179 return atoms
1182# In the interest of using the same RE for both atomic positions and forces,
1183# we make one of the columns optional. That's because atomic positions have
1184# 6 columns, while forces only has 5 columns. Otherwise they are very similar.
1185_re_atom = re.compile(
1186 r'^\s*\S+\s+(\S+)\s+(?:\S+\s+)?(\S+)\s+(\S+)\s+(\S+)\s*$'
1187)
1188_re_forceblock = re.compile(r'^\s*Center\s+Atomic\s+Forces\s+\S+\s*$')
1189_re_l716 = re.compile(r'^\s*\(Enter .+l716.exe\)$')
1192def _compare_merge_configs(configs, new):
1193 """Append new to configs if it contains a new geometry or new data.
1195 Gaussian sometimes repeats a geometry, for example at the end of an
1196 optimization, or when a user requests vibrational frequency
1197 analysis in the same calculation as a geometry optimization.
1199 In those cases, rather than repeating the structure in the list of
1200 returned structures, try to merge results if doing so doesn't change
1201 any previously calculated values. If that's not possible, then create
1202 a new "image" with the new results.
1203 """
1204 if not configs:
1205 configs.append(new)
1206 return
1208 old = configs[-1]
1210 if old != new:
1211 configs.append(new)
1212 return
1214 oldres = old.calc.results
1215 newres = new.calc.results
1216 common_keys = set(oldres).intersection(newres)
1218 for key in common_keys:
1219 if np.any(oldres[key] != newres[key]):
1220 configs.append(new)
1221 return
1222 else:
1223 oldres.update(newres)
1226def read_gaussian_out(fd, index=-1):
1227 configs = []
1228 atoms = None
1229 energy = None
1230 dipole = None
1231 forces = None
1232 for line in fd:
1233 line = line.strip()
1234 if line.startswith(r'1\1\GINC'):
1235 # We've reached the "archive" block at the bottom, stop parsing
1236 break
1238 if (line == 'Input orientation:'
1239 or line == 'Z-Matrix orientation:'):
1240 if atoms is not None:
1241 atoms.calc = SinglePointCalculator(
1242 atoms, energy=energy, dipole=dipole, forces=forces,
1243 )
1244 _compare_merge_configs(configs, atoms)
1245 atoms = None
1246 energy = None
1247 dipole = None
1248 forces = None
1250 numbers = []
1251 positions = []
1252 pbc = np.zeros(3, dtype=bool)
1253 cell = np.zeros((3, 3))
1254 npbc = 0
1255 # skip 4 irrelevant lines
1256 for _ in range(4):
1257 fd.readline()
1258 while True:
1259 match = _re_atom.match(fd.readline())
1260 if match is None:
1261 break
1262 number = int(match.group(1))
1263 pos = list(map(float, match.group(2, 3, 4)))
1264 if number == -2:
1265 pbc[npbc] = True
1266 cell[npbc] = pos
1267 npbc += 1
1268 else:
1269 numbers.append(max(number, 0))
1270 positions.append(pos)
1271 atoms = Atoms(numbers, positions, pbc=pbc, cell=cell)
1272 elif (line.startswith('Energy=')
1273 or line.startswith('SCF Done:')):
1274 # Some semi-empirical methods (Huckel, MINDO3, etc.),
1275 # or SCF methods (HF, DFT, etc.)
1276 energy = float(line.split('=')[1].split()[0].replace('D', 'e'))
1277 energy *= Hartree
1278 elif (line.startswith('E2 =') or line.startswith('E3 =')
1279 or line.startswith('E4(') or line.startswith('DEMP5 =')
1280 or line.startswith('E2(')):
1281 # MP{2,3,4,5} energy
1282 # also some double hybrid calculations, like B2PLYP
1283 energy = float(line.split('=')[-1].strip().replace('D', 'e'))
1284 energy *= Hartree
1285 elif line.startswith('Wavefunction amplitudes converged. E(Corr)'):
1286 # "correlated method" energy, e.g. CCSD
1287 energy = float(line.split('=')[-1].strip().replace('D', 'e'))
1288 energy *= Hartree
1289 elif _re_l716.match(line):
1290 # Sometimes Gaussian will print "Rotating derivatives to
1291 # standard orientation" after the matched line (which looks like
1292 # "(Enter /opt/gaussian/g16/l716.exe)", though the exact path
1293 # depends on where Gaussian is installed). We *skip* the dipole
1294 # in this case, because it might be rotated relative to the input
1295 # orientation (and also it is numerically different even if the
1296 # standard orientation is the same as the input orientation).
1297 line = fd.readline().strip()
1298 if not line.startswith('Dipole'):
1299 continue
1300 dip = line.split('=')[1].replace('D', 'e')
1301 tokens = dip.split()
1302 dipole = []
1303 # dipole elements can run together, depending on what method was
1304 # used to calculate them. First see if there is a space between
1305 # values.
1306 if len(tokens) == 3:
1307 dipole = list(map(float, tokens))
1308 elif len(dip) % 3 == 0:
1309 # next, check if the number of tokens is divisible by 3
1310 nchars = len(dip) // 3
1311 for i in range(3):
1312 dipole.append(float(dip[nchars * i:nchars * (i + 1)]))
1313 else:
1314 # otherwise, just give up on trying to parse it.
1315 dipole = None
1316 continue
1317 # this dipole moment is printed in atomic units, e-Bohr
1318 # ASE uses e-Angstrom for dipole moments.
1319 dipole = np.array(dipole) * Bohr
1320 elif _re_forceblock.match(line):
1321 # skip 2 irrelevant lines
1322 fd.readline()
1323 fd.readline()
1324 forces = []
1325 while True:
1326 match = _re_atom.match(fd.readline())
1327 if match is None:
1328 break
1329 forces.append(list(map(float, match.group(2, 3, 4))))
1330 forces = np.array(forces) * Hartree / Bohr
1331 if atoms is not None:
1332 atoms.calc = SinglePointCalculator(
1333 atoms, energy=energy, dipole=dipole, forces=forces,
1334 )
1335 _compare_merge_configs(configs, atoms)
1336 return configs[index]