Gaussian

Gaussian is a computational chemistry code based on gaussian basis functions.

Setup

Ask your system administrator to install Gaussian for you.

The ASE Gaussian calculator has been written with Gaussian 16 (g16) in mind, but it will likely work with newer and older versions of Gaussian as well. By default, the Calculator will look for executables named g16, g09, and g03 in that order. If your Gaussian executable is named differently, or if it is not present in PATH, then you must pass the path and name of your Gaussian executable to the command keyword argument of the Gaussian calculator. The default command looks like g16 < PREFIX.com > PREFIX.log, so template the command similarly. Alternatively, you may set the ASE_GAUSSIAN_COMMAND environment variable to the full Gaussian executable command.

Examples

Here is a command line example of how to optimize the geometry of a water molecule using the PBE density functional:

$ ase build H2O | ase run gaussian -p xc=PBE,basis=3-21G -f 0.02
$ ase gui stdin.traj@-1 -tg "a(1,0,2),d(0,1)"
102.58928991353669 1.0079430292939233

An example of creating a Gaussian calculator in the python interface is:

from ase.calculators.gaussian import Gaussian

calc = Gaussian(label='calc/gaussian',
                xc='B3LYP',
                basis='6-31+G*',
                scf='maxcycle=100')

Parameters

The Gaussian calculator has three main types of parameters:

  1. Link0 keywords

  2. Route section keywords

  3. ASE-specific keywords, or convenience keywords.

The Gaussian calculator maintains a list of Link0 keywords and ASE-specific keywords. Any keyword not on one of those two lists is assumed to be a route section keyword, and will be placed in the Gaussian input file accordingly.

For example, consider the following Gaussian input file:

%mem=1GB
%chk=MyJob.chk
%save
#P b3lyp/6-31G scf=qc

My job label

0 1
H 0.00 0.00 0.00
H 0.00 0.00 0.74

This would be generated with the following Python code:

from ase import Atoms
from ase.calculators.gaussian import Gaussian

atoms = Atoms('H2', [[0, 0, 0], [0, 0, 0.74]])
atoms.calc = Gaussian(mem='1GB',
                      chk='MyJob.chk',
                      save=None,
                      method='b3lyp',
                      basis='6-31G',
                      scf='qc')
atoms.get_potential_energy()

Alternatively, you may use the xc keyword in place of the method keyword. xc is almost identical to method, except that xc can translate between the common definitions of some exchange-correlation functionals and Gaussian’s name for those functions, for example PBE to PBEPBE. The method keyword will not do any translation, whatever value you provide to method will be written to the input file verbatim. If both are provided, method overrides xc.

Note that the Gaussian calculator puts each route keyword on its own line, though this should not affect the result of the calculation.

When a route section keyword has multiple arguments, it is usually written like scf(qc,maxcycle=1000) in the Gaussian input file. There are at least two ways of generating this with the Gaussian calculator: Gaussian(scf="qc,maxcycle=100") and Gaussian(scf=['qc', 'maxcycle=100']), with the latter being somewhat more convenient for scripting purposes.

Aside from the link-line and route section arguments, the Gaussian calculator accepts a few additional convenience arguments.

keyword

type

default value

description

label

str

'Gaussian'

Name to use for input and output files.

output_type

str

P

Level of output to record in the Gaussian output file - this may be N- normal or P - additional.

method

str

None

Level of theory to use, e.g. hf, ccsd, mp2, or b3lyp. Overrides xc (see below).

xc

str

None

Level of theory to use. Translates several XC functionals from their common name (e.g. PBE) to their internal Gaussian name (e.g. PBEPBE).

basis

str

None

The basis set to use. If not provided, no basis set will be requested, which usually results in STO-3G. Maybe omitted if basisfile is set (see below).

fitting_basis

str

None

The fitting basis set to use.

charge

int

See description

The system charge. If not provided, it will be automatically determined from the Atoms object’s initial_charges.

mult

int

See description

The system multiplicity (spin + 1). If not provided, it will be automatically determined from the Atoms object’s initial_magnetic_moments.

basisfile

str

None

The basis file to use. If a value is provided, basis may be omitted (it will be automatically set to 'gen')

basis_set

str

None

The basis set definition to use. This is an alternative to basisfile, and would be the same as the contents of such a file.

extra

str

None

Extra lines to be included in the route section verbatim. It should not be necessary to use this, but it is included for backwards compatibility.

addsec

str

None

Text to be added after the molecular geometry specification, e.g. for defining constraints with opt='modredundant'.

ioplist

list

None

A collection of IOPs definitions to be included in the route line.

spinlist

list

None

A list of nuclear spins to be added into the nuclear properties section of the molecule specification.

zefflist

list

None

A list of effective charges to be added into the nuclear properties section of the molecule specification.

qmomlist

list

None

A list of nuclear quadropole moments to be added into the nuclear properties section of the molecule specification.

nmagmlist

list

None

A list of nuclear magnetic moments to be added into the nuclear properties section of the molecule specification.

znuclist

list

None

A list of nuclear charges to be added into the nuclear properties section of the molecule specification.

radnuclearlist

list

None

A list of nuclear radii to be added into the nuclear properties section of the molecule specification.

GaussianOptimizer and GaussianIRC

There are also two Gaussian-specific Optimizer-like classes: GaussianOptimizer and GaussianIRC, which can be used for geometry optimizations and IRC calculations, respectively. These can be invoked in the following way:

from ase.calculators.gaussian import Gaussian, GaussianOptimizer
atoms = ...
calc_opt = Gaussian(...)
opt = GaussianOptimizer(atoms, calc_opt)
opt.run(fmax='tight', steps=100)

Note that this differs from ASE’s standard Optimizer classes in a few key ways:

  1. The fmax keyword takes a string rather than a force/energy criterion. Valid keywords are described in the Gaussian manual page for optimization.

  2. Unlike ASE’s standard Optimizer classes, it is not possible to iterate over the optimization with opt.irun(...).

  3. It is also not possible to create a Trajectory file which records the optimization with opt = GaussianOptimizer(..., trajectory='opt.traj'). However, it should be possible to obtain the trajectory by reading the Gaussian output file after the optimization has finished.

Additional arguments to Gaussian’s opt keyword can be passed to the calculator in the following way:

opt.run(fmax='tight', steps=100, opt='calcfc,ts')

This example requests a Hessian calculation followed by optimization to a saddle point (“transition state optimization”).

The GaussianIRC class can also be used to run IRC or pseudo-IRC calculations. For example, the following script optimizes to a saddle point, then runs an IRC optimization in the forward- and reverse-direction:

from ase.calculators.gaussian import Gaussian, GaussianOptimizer, GaussianIRC
atoms = ...

# Optimize to a saddle point
calc_opt = Gaussian(label='opt', ...)
opt = GaussianOptimizer(atoms, calc_opt)
opt.run(fmax='tight', steps=100, opt='calcfc,ts')
tspos = atoms.positions.copy()

# Do a vibrational frequency calculation and store the Hessian in a
# checkpoint file, for use in subsequent IRC calculations
atoms.calc = Gaussian(label='sp', chk='sp.chk', freq='')
atoms.get_potential_energy()

# Perform IRC in the "forwards" direction
calc_irc_for = Gaussian(label='irc_for', chk='irc_for.chk', oldchk='sp.chk', ...)
irc_for = GaussianIRC(atoms, calc_irc_for)
irc_for.run(direction='forward', steps=20, irc='rcfc')  # reuses Hessian

# Perform IRC in the "reverse" direction
# First, restore TS positions
atoms.positions[:] = tspos
calc_irc_rev = Gaussian(label='irc_rev', chk='irc_rev.chk', oldchk='sp.chk', ...)
irc_rev = GaussianIRC(atoms, calc_irc_rev)
irc_rev.run(direction='reverse', steps=20, irc='rcfc')

It should also be possible to use the same Gaussian calculator object for each of these steps, so long as the label is changed between calculations (to avoid overwriting the output file) and the settings are changed appropriately. It should also be possible to use the same GaussianIRC object for both the forwards and reverse IRC calculations, so long as the label is changed (again to avoid overwriting the output file).

class ase.calculators.gaussian.Gaussian(*args, label='Gaussian', **kwargs)[source]

File-IO calculator.

command: str

Command used to start calculation.