6. Calculators

Exploring potential energy surfaces (PESs) requires calculation of energies and its derivatives (gradients, Hessian matrices).

For testing purposes and method development, pysisyphus implements various analytical 2d potentials, as they allow fast evaluations of the aforementioned quantities. Actual (production) calculations are carried out by wrapping existing quantum chemistry codes (ORCA, TURBOMOLE, Gaussian etc.) and delegating the required calculations to them. Pysisyphus generates all necessary inputs, executes the QC code and parses their output for the requested quantities.

Furthermore pysisyphus provides several "meta"-calculators which wrap other (actual) calculators, to modify calculated energies and forces. Examples for this are the Dimer calculator, used for carrying out transition state searches or the ONIOM calculator, allowing multi-level calculations comprising different levels of theory.

External forces, e.g. a restraining spherical potential or harmonic restraints on primitive internal coordinates (stretches, bends, torsion) can be applied with ExternalPotential.

6.1. YAML input

Possible keywords for the YAML input can be derived from inspecting the possible arguments of the Calculator base-class (see below) and the possible arguments of the respective calculator, e.g. ORCA or XTB.

The most commonly used keywords, derived from the Calculator baseclass are mem, handling the requested memory per core in MB, pal, handling the number of requested CPU cores, charge, the total charge of the system and mult, the systems multiplicity.

For (excited state, ES) calculations carried out with calculators derived from OverlapCalculator additional keywords are possible. The most common keywords controlling ES calculations are track, activating ES tracking, ovlp_type, selecting the tracking algorithm and ovlp_with, handling the selection of the reference cycle.

An example input highlighting the most important keywords is shown below.

geom:
 [... omitted ...]
calc:
 type: orca                     # Calculator type, e.g. g09/g16/openmolcas/
                                # orca/pyscf/turbomole/dftb+/mopac/psi4/xtb

 pal: 1                         # Number of CPU cores
 mem: 1000                      # Memory per core
 charge: 0                      # Charge
 mult: 1                        # Multiplicity
 # Keywords for ES calculations
 track: False                   # Activate ES tracking
 ovlp_type: tden                # Tracking algorithm
 ovlp_with: previous            # Reference cycle selection
 # Additional calculator specific keywords
 [... omitted ...]
opt:
 [... omitted ...]

6.2. Calculator base classes

class pysisyphus.calculators.Calculator.Calculator(calc_number=0, charge=0, mult=1, base_name='calculator', pal=1, mem=1000, check_mem=True, retry_calc=1, last_calc_cycle=None, clean_after=True, out_dir='qm_calcs')

__init__(calc_number=0, charge=0, mult=1, base_name='calculator', pal=1, mem=1000, check_mem=True, retry_calc=1, last_calc_cycle=None, clean_after=True, out_dir='qm_calcs')

Base-class of all calculators.

Meant to be extended.

Parameters

• calc_number (int, default=0) -- Identifier of the Calculator. Used in distinguishing it from other Calculators, e.g. in ChainOfStates calculations. Also used in the creation of filenames.

• charge (int, default=0) -- Molecular charge.

• mult (int, default=1) -- Molecular multiplicity (1 = singlet, 2 = doublet, ...)

• base_name (str, default=calculator) -- Generated filenames will start with this string.

• pal (int, default=1) -- Positive integer that gives the number of physical cores to use on 1 node.

• mem (int, default=1000) -- Mememory per core in MB. The total amount of memory is given as mem*pal.

• check_mem (bool, default=True) -- Whether to adjust the requested memory if too much is requested.

• retry_calc (int, default=0) -- Number of additional retries when calculation failed.

• last_calc_cycle (int) -- Internal variable used in restarts.

• clean_after (bool) -- Delete temporary directory were calculations were executed after a calculation.

• out_dir (str) -- Path that is prepended to generated filenames.

clean(path)

Delete the temporary directory.

Parameters

path (Path) -- Directory to delete.

conf_key = None

get_cmd(key='cmd')

get_energy(atoms, coords)

Meant to be extended.

get_forces(atoms, coords)

Meant to be extended.

get_hessian(atoms, coords)

Meant to be extended.

get_restart_info()

Return a dict containing chkfiles.

Returns

restart_info -- Dictionary holding the calculator state. Used for restoring calculaters in restarted calculations.

Return type

dict

keep(path)

Backup calculation results.

Parameters

path (Path) -- Temporary directory of the calculation.

Returns

kept_fns -- Dictonary holding the filenames that were backed up. The keys correspond to the type of file.

Return type

dict

log(message='')

Write a log message.

Wraps the logger variable.

Parameters

message (str) -- Message to be logged.

make_fn(name, counter=None, return_str=False)

Make a full filename.

Return a full filename including the calculator name and the current counter given a suffix.

Parameters

• name (str) -- Suffix of the filename.

• counter (int, optional) -- If not given use the current calc_counter.

• return_str (int, optional) -- Return a string instead of a Path when True.

Returns

fn -- Filename.

Return type

str

property name

prepare(inp)

Prepare a temporary directory and write input.

Similar to prepare_path, but the input is also written into the prepared directory.

inpstr

Input to be written into the file self.inp_fn in the prepared directory.

Returns

path -- Prepared directory.

Return type

Path

prepare_coords(atoms, coords)

Get 3d coords in Angstrom.

Reshape internal 1d coords to 3d and convert to Angstrom.

Parameters

• atoms (iterable) -- Atom descriptors (element symbols).

• coords (np.array, 1d) -- 1D-array holding coordinates in Bohr.

Returns

coords -- 3D-array holding coordinates in Angstrom.

Return type

np.array, 3d

prepare_input(atoms, coords, calc_type)

Meant to be extended.

prepare_path(use_in_run=False)

Get a temporary directory handle.

Create a temporary directory that can later be used in a calculation.

Parameters

use_in_run (bool, option) -- Sets the internal variable self.path_already_prepared that is later read by self.run(). No new temporary directory will be created in self.run().

Returns

path -- Prepared directory.

Return type

Path

prepare_pattern(raw_pat)

Prepare globs.

Transforms an entry of self.to_keep into a glob and a key suitable for the use in self.keep().

Parameters

raw_pat (str) -- Entry of self.to_keep

Returns

• pattern (str) -- Glob that can be used in Path.glob()

• multi (bool) -- Flag if glob may match multiple files.

• key (str) -- A key to be used in the kept_fns dict.

prepare_turbo_coords(atoms, coords)

Get a Turbomole coords string.

Parameters

• atoms (iterable) -- Atom descriptors (element symbols).

• coords (np.array, 1d) -- 1D-array holding coordinates in Bohr.

Returns

coords -- String holding coordinates in Turbomole coords format.

Return type

str

prepare_xyz_string(atoms, coords)

Returns a xyz string in Angstrom.

Parameters

• atoms (iterable) -- Atom descriptors (element symbols).

• coords (np.array, 1d) -- 1D-array holding coordinates in Bohr.

Returns

xyz_str -- Coordinates in .xyz format.

Return type

string

print_capabilities()

print_out_fn(path)

Print calculation output.

Prints the output of a calculator after a calculation.

Parameters

path (Path) -- Temporary directory of the calculation.

run(inp, calc, add_args=None, env=None, shell=False, hold=False, keep=True, cmd=None, inc_counter=True, run_after=True, parser_kwargs=None, symlink=True)

Run a calculation.

The bread-and-butter method to actually run an external quantum chemistry code.

Parameters

• inp (str) -- Input for the external program that is written to the temp-dir.

• calc (str, hashable) -- Key (and more or less type of calculation) to select the right parsing function from self.parser_funcs.

• add_args (iterable, optional) -- Additional arguments that will be appended to the program call.

• env (Environment, optional) -- A potentially modified environment for the subprocess call.

• shell (bool, optional) -- Use a shell to execute the program call. Need for Turbomole were we chain program calls like dscf; escf.

• hold (bool, optional) -- Wether to remove the temporary directory after the calculation.

• keep (bool, optional) -- Wether to backup files as specified in self.to_keep(). Usually you want this.

• cmd (str or iterable, optional) -- Overwrites self.base_cmd.

• inc_counter (bool, optional) -- Wether to increment the counter after a calculation.

Returns

results -- Dictionary holding all applicable results of the calculations like the energy, a forces vector and/or excited state energies from TDDFT.

Return type

dict

run_after(path)

Meant to be extended.

This method is called after a calculation was done, but before entering self.keep() and self.clean(). Can be used to call tools like formchk or ricctools.

set_restart_info(restart_info)

Sets restart information (chkfiles etc.) on the calculator.

Parameters

restart_info (dict) -- Dictionary holding the calculator state. Used for restoring calculaters in restarted calculations.

verify_chkfiles(chkfiles)

Checks if given chkfiles exist and return them as Paths

Parameters

chkfiles (dict) -- Dictionary holding the chkfiles. The keys correspond to the attribute names, the values are strs holding the (potentially full) filename (path).

Returns

paths -- Dictionary of Paths.

Return type

dict

6.3. OverlapCalculator base class

class pysisyphus.calculators.OverlapCalculator.NTOs(ntos, lambdas)

Bases: tuple

lambdas

Alias for field number 1

ntos

Alias for field number 0

class pysisyphus.calculators.OverlapCalculator.OverlapCalculator(*args, track=False, ovlp_type='tden', double_mol=False, ovlp_with='previous', XY='X+Y', adapt_args=(0.5, 0.3, 0.6), use_ntos=4, pr_nto=False, nto_thresh=0.3, cdds=None, orient='', dump_fn='overlap_data.h5', h5_dump=False, ncore=0, conf_thresh=0.001, dyn_roots=0, mos_ref='cur', mos_renorm=True, **kwargs)

Bases: Calculator

H5_MAP = {'all_energies': 'all_energies_list', 'ci_coeffs': 'ci_coeff_list', 'coords': 'coords_list', 'mo_coeffs': 'mo_coeff_list', 'ref_roots': 'reference_roots', 'roots': 'roots_list'}

OVLP_TYPE_VERBOSE = {'nto': 'natural transition orbital overlap', 'nto_org': 'original natural transition orbital overlap', 'tden': 'transition density matrix overlap', 'wf': 'wavefunction overlap'}

VALID_CDDS = (None, 'calc', 'render')

VALID_KEYS = ['wf', 'tden', 'nto', 'nto_org']

VALID_XY = ('X', 'X+Y', 'X-Y')

blowup_ci_coeffs(ci_coeffs)

calc_cdd_cube(root, cycle=- 1)

calculate_state_ntos(state_ci_coeffs, mos)

property data_model

dump_overlap_data()

static from_overlap_data(h5_fn, set_wfow=False)

get_h5_group()

get_indices(indices=None)

A new root is determined by selecting the overlap matrix row corresponding to the reference root and checking for the root with the highest overlap (at the current geometry).

The overlap matrix is usually formed by a double loop like:

overlap_matrix = np.empty((ref_states, cur_states)) for i, ref_state in enumerate(ref_states):

for j, cur_state in enumerate(cur_states):

overlap_matrix[i, j] = make_overlap(ref_state, cur_state)

So the reference states run along the rows. Thats why the ref_state index comes first in the 'indices' tuple.

static get_mo_norms(mo_coeffs, ao_ovlp)

MOs are in rows.

get_nto_overlaps(indices=None, ao_ovlp=None, org=False)

get_orbital_matrices(indices=None, ao_ovlp=None)

Return MO coefficents and AO overlaps for the given indices.

If not provided, a AO overlap matrix is constructed from one of the MO coefficient matrices (controlled by self.mos_ref). Also, if requested one of the two MO coefficient matrices is re-normalized.

get_ref_mos(ref_mo_coeffs, cur_mo_coeffs)

static get_sao_from_mo_coeffs(mo_coeffs)

Recover AO overlaps from given MO coefficients.

For MOs in the columns of mo_coeffs:

S_AO = C⁻¹^T C⁻¹ S_AO C = C⁻¹^T (S_AO C)^T = C⁻¹ C^T S_AO^T = C⁻¹ C^T S_AO C = I

Here, MOs are expected to be in rows of mo_coeffs, yielding

C S_AO C^T = I

get_sao_from_mo_coeffs_and_dump(mo_coeffs)

get_tden_overlaps(indices=None, ao_ovlp=None)

get_wf_overlaps(indices=None, ao_ovlp=None)

nto_org_overlaps(ntos_1, ntos_2, ao_ovlp, nto_thresh=0.3)

nto_overlaps(ntos_1, ntos_2, ao_ovlp)

prepare_overlap_data(path)

Implement calculator specific parsing of MO coefficients and CI coefficients here. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

render_cdd_cube()

static renorm_mos(mo_coeffs, ao_ovlp)

property roots_number

set_ntos(mo_coeffs, ci_coeffs)

set_wfow(ci_coeffs)

store_overlap_data(atoms, coords, path=None, overlap_data=None)

tden_overlap_with_calculator(calc, ao_ovlp=None)

tden_overlaps(mo_coeffs1, ci_coeffs1, mo_coeffs2, ci_coeffs2, ao_ovlp)

Parameters

• mo_coeffs1 (ndarray, shape (MOs, AOs)) -- MO coefficient matrix. One row per MO, one column per basis function. Usually square.

• mo_coeffs2 (ndarray) -- See mo_coeffs1.

• ci_coeffs1 (ndarray, shape(occ. MOs, MOs)) -- CI-coefficient matrix.

• ci_coeffs2 (ndarray, shape(occ. MOs, MOs)) -- See ci_coeffs1.

• ao_ovlp (ndarray, shape(AOs1, AOs2)) -- Double molcule AO overlaps.

track_root(ovlp_type=None)

Check if a root flip occured occured compared to the previous cycle by calculating the overlap matrix wrt. a reference cycle.

update_array_dims()

Assure consistent array shapes when the number of calculated roots changed.

wf_overlap_with_calculator(calc, ao_ovlp=None)

wf_overlaps(mo_coeffs1, ci_coeffs1, mo_coeffs2, ci_coeffs2, ao_ovlp=None)

pysisyphus.calculators.OverlapCalculator.get_data_model(exc_state_num, occ_mo_num, virt_mo_num, ovlp_type, atoms, max_cycles)

6.4. Calculators with Excited state capabilities

6.4.1. Gaussian09

class pysisyphus.calculators.Gaussian09.Gaussian09(*args, **kwargs)

conf_key = 'gaussian09'

6.4.2. Gaussian16

class pysisyphus.calculators.Gaussian16.Gaussian16(route, gbs='', gen='', keep_chk=False, stable='', fchk=None, **kwargs)

Bases: OverlapCalculator

conf_key = 'gaussian16'

get_chkfiles()

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)

Meant to be extended.

keep(path)

Backup calculation results.

Parameters

path (Path) -- Temporary directory of the calculation.

Returns

kept_fns -- Dictonary holding the filenames that were backed up. The keys correspond to the type of file.

Return type

dict

make_exc_str()

make_fchk(path)

make_gbs_str()

parse_635r_dump(dump_path, roots, nmos)

parse_charges(path=None)

parse_double_mol(path, out_fn=None)

static parse_fchk(fchk_path, keys)

parse_force(path)

parse_hessian(path)

parse_keyword(text)

parse_log(*args, **kwargs)

parse_stable(path)

parse_tddft(path)

prepare_input(atoms, coords, calc_type, did_stable=False, point_charges=None)

Meant to be extended.

prepare_overlap_data(path)

Implement calculator specific parsing of MO coefficients and CI coefficients here. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

reuse_data(path)

run_after(path)

Meant to be extended.

This method is called after a calculation was done, but before entering self.keep() and self.clean(). Can be used to call tools like formchk or ricctools.

run_calculation(atoms, coords, **prepare_kwargs)

run_double_mol_calculation(atoms, coords1, coords2)

run_rwfdump(path, rwf_index, chk_path=None)

run_stable(atoms, coords, **prepare_kwargs)

set_chkfiles(chkfiles)

6.4.3. OpenMolcas

class pysisyphus.calculators.OpenMolcas.OpenMolcas(basis, inporb, roots, mdrlxroot, supsym=None, track=True, **kwargs)

Bases: Calculator

build_rassi_str()

build_supsym_str(supsym)

Can handle only one subgroup for now.

conf_key = 'openmolcas'

get_energy(atoms, coords)

Meant to be extended.

get_forces(atoms, coords)

Meant to be extended.

get_pal_env()

keep(path)

Backup calculation results.

Parameters

path (Path) -- Temporary directory of the calculation.

Returns

kept_fns -- Dictonary holding the filenames that were backed up. The keys correspond to the type of file.

Return type

dict

parse_energies(text)

parse_gradient(path)

parse_rassi_track(path)

prepare_coords(atoms, coords)

Get 3d coords in Angstrom.

Reshape internal 1d coords to 3d and convert to Angstrom.

Parameters

• atoms (iterable) -- Atom descriptors (element symbols).

• coords (np.array, 1d) -- 1D-array holding coordinates in Bohr.

Returns

coords -- 3D-array holding coordinates in Angstrom.

Return type

np.array, 3d

prepare_input(atoms, coords)

Meant to be extended.

reattach(last_calc_cycle)

6.4.4. ORCA 4.2.1 / 5.0.1

class pysisyphus.calculators.ORCA.ORCA(keywords, blocks='', gbw=None, do_stable=False, numfreq=False, **kwargs)

Bases: OverlapCalculator

__init__(keywords, blocks='', gbw=None, do_stable=False, numfreq=False, **kwargs)

ORCA calculator.

Wrapper for creating ORCA input files for energy, gradient and Hessian calculations. The PAL and memory inputs must not be given in the keywords and/or blocks, as they are handled by the 'pal' and 'memory' arguments.

Parameters

• keywords (str) -- Keyword line, as normally given in ORCA, excluding the leading "!".

• blocks (str, optional) -- ORCA block input(s), e.g. for TD-DFT calculations (%tddft ... end). As the blocks start with a leading "%", wrapping the input in quotes ("") is required, otherwise the parsing will fail.

• gbw (str, optional) -- Path to an input gbw file, which will be used as initial guess for the first calculation. Will be overriden later, with the path to the gbw file of a previous calculation.

• do_stable (bool, optional) -- Run stability analysis until a stable wavefunction is obtained, before every calculation.

• numfreq (boo, optional) -- Use numerical frequencies instead of analytical ones.

• mem (int) -- Mememory per core in MB.

check_termination(*args, **kwargs)

clean_tmp(path)

conf_key = 'orca'

get_block_str()

get_chkfiles()

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_moinp_str(gbw)

get_stable_wavefunction(atoms, coords)

keep(path)

Backup calculation results.

Parameters

path (Path) -- Temporary directory of the calculation.

Returns

kept_fns -- Dictonary holding the filenames that were backed up. The keys correspond to the type of file.

Return type

dict

parse_all_energies(text=None, triplets=None)

static parse_atoms_coords(inp, *args, **kwargs)

parse_cis(cis)

Read binary CI vector file from ORCA.

Adapted from TheoDORE 1.7.1, Authors: S. Mai, F. Plasser https://sourceforge.net/p/theodore-qc

parse_energy(path)

parse_engrad(path)

static parse_engrad_info(inp, *args, **kwargs)

parse_gbw(gbw_fn)

Adapted from https://orcaforum.kofo.mpg.de/viewtopic.php?f=8&t=3299

The first 5 long int values represent pointers into the file:

Pointer @+0: Internal ORCA data structures Pointer @+8: Geometry Pointer @+16: BasisSet Pointer @+24: Orbitals Pointer @+32: ECP data

static parse_hess_file(inp, *args, **kwargs)

parse_hessian(path)

parse_mo_numbers(out_fn)

parse_stable(path)

prepare_input(atoms, coords, calc_type, point_charges=None, do_stable=False)

Meant to be extended.

prepare_overlap_data(path)

Implement calculator specific parsing of MO coefficients and CI coefficients here. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

reattach(last_calc_cycle)

run_after(path)

Meant to be extended.

This method is called after a calculation was done, but before entering self.keep() and self.clean(). Can be used to call tools like formchk or ricctools.

run_calculation(atoms, coords, **prepare_kwargs)

Basically some kind of dummy method that can be called to execute ORCA with the stored cmd of this calculator.

set_chkfiles(chkfiles)

set_mo_coeffs(mo_coeffs=None, gbw=None)

static set_mo_coeffs_in_gbw(in_gbw_fn, out_gbw_fn, mo_coeffs)

See self.parse_gbw.

store_and_track(results, func, atoms, coords, **prepare_kwargs)

pysisyphus.calculators.ORCA.make_sym_mat(table_block)

pysisyphus.calculators.ORCA.save_orca_pc_file(point_charges, pc_fn, hardness=None)

6.4.5. PySCF 1.7.6

6.4.6. Turbomole 7.x

Pysisyphus does not implement a wrapper for define, so the user has to manually prepare a directory containing a valid control file. An automated define wrapper, restricted to ground state functionality, is available via the QCEngine project, to which I contributed the Turbomole harness.

Care should be taken to include only the minimum amount of necessary files in the control_path directory, e.g., (auxbasis, basis, control, coord, mos) for a closed-shell calculation using RI. A gradient file must not be present in control_path, as well as other subdirectories and files with .out extension. The coord file, while not strictly required, should be kept too, to facilitate testing of the setup with standalone Turbomole.

It may be a good idea to pre-converge the calculation in control_path, to see if the setup is correct and actually works. Resulting files like energy, statistics can be deleted; mos should be kept, as the converged MOs are reused in pysisyphus.

If an excited-state optimization is desired, care has to be taken, to include $exopt [n] for TD-DFT/TDA or the geoopt state=([n]) (ricc2)! Tracking of excited states is currently possible for closed shell egrad and ricc2 calculations.

The current implementation was tested against Turbomole 7.4.1 and QCEngine 0.19.0. Please see examples/complex/11_turbomole_gs_tsopt for a full example where Turbmole is utilized in a growing string calculation. The same example, using QCEngine, is found in examples/complex/12_qcengine_turbomole_gs_tsopt. MOs are not reused with the QCEngine calculator, so the native pysisyphus calculator is probably faster.

class pysisyphus.calculators.Turbomole.Turbomole(control_path, root=None, double_mol_path=None, **kwargs)

Bases: OverlapCalculator

append_control(to_append, log_msg='', **kwargs)

conf_key = 'turbomole'

get_chkfiles()

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_forces(atoms, coords, cmd=None, **prepare_kwargs)

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_pal_env()

get_ricc2_root(text)

keep(path)

Backup calculation results.

Parameters

path (Path) -- Temporary directory of the calculation.

Returns

kept_fns -- Dictonary holding the filenames that were backed up. The keys correspond to the type of file.

Return type

dict

make_molden(path)

parse_cc2_vectors(ccre)

parse_double_mol(path)

Parse a double molecule overlap matrix from Turbomole output to be used with WFOWrapper.

parse_energy(path)

parse_force(path)

parse_gs_energy()

Several places are possible: $subenergy from control file total energy from turbomole.out Final MP2 energy from turbomole.out with ADC(2) Final CC2 energy from turbomole.out with CC(2)

parse_hessian(path, fn=None)

parse_mos()

parse_td_vectors(text)

For TDA calculations only the X vector is present in the ciss_a/etc. file. In TDDFT calculations there are twise as much items compared with TDA. The first half corresponds to (X+Y) and the second half to (X-Y). X can be calculated as X = ((X+Y)+(X-Y))/2. Y is then given as Y = (X+Y)-X. The normalization can then by checked as

np.concatenate((X, Y)).dot(np.concatenate((X, -Y)))

and should be 1.

prepare_input(atoms, coords, calc_type, point_charges=None)

To rectify this we have to construct the basecmd dynamically and construct it ad hoc. We could set a RI flag in the beginning and select the correct scf binary here from it. Then we select the following binary on demand, e.g. aoforce or rdgrad or egrad etc.

prepare_overlap_data(path)

Implement calculator specific parsing of MO coefficients and CI coefficients here. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

prepare_point_charges(point_charges)

$point_charges <x> <y> <z> <q>

prepare_td(text)

run_after(path)

Meant to be extended.

This method is called after a calculation was done, but before entering self.keep() and self.clean(). Can be used to call tools like formchk or ricctools.

run_calculation(atoms, coords)

Basically some kind of dummy method that can be called to execute Turbomole with the stored cmd of this calculator.

run_double_mol_calculation(atoms, coords1, coords2)

set_chkfiles(chkfiles)

set_occ_and_mo_nums(text)

store_and_track(results, func, atoms, coords, **prepare_kwargs)

sub_control(pattern, repl, log_msg='', **kwargs)

6.4.7. DFTB+ 20.x

class pysisyphus.calculators.DFTBp.DFTBp(parameter, *args, slakos=None, root=None, **kwargs)

Bases: OverlapCalculator

conf_key = 'dftbp'

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

static get_excited_state_str(root, forces=False)

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

static get_gen_str(atoms, coords)

max_ang_moms = {'mio-ext': {'C': 'p', 'H': 's', 'N': 'p', 'O': 'p'}}

parse_energy(path)

static parse_exc_dat(text)

parse_forces(path)

parse_total_energy(text)

prepare_input(atoms, coords, calc_type)

Meant to be extended.

prepare_overlap_data(path)

Implement calculator specific parsing of MO coefficients and CI coefficients here. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

run_calculation(atoms, coords, **prepare_kwargs)

pysisyphus.calculators.DFTBp.parse_mo(eigvec)

pysisyphus.calculators.DFTBp.parse_xplusy(text)

6.5. Calculators with Ground state capabilities

6.5.1. MOPAC 2016

class pysisyphus.calculators.MOPAC.MOPAC(method='PM7', **kwargs)

Bases: Calculator

http://openmopac.net/manual/

CALC_TYPES = {'energy': '1SCF', 'gradient': '1SCF GRADIENTS', 'hessian': 'DFORCE FORCE LET'}

METHODS = ['am1', 'pm3', 'pm6', 'pm6-dh2', 'pm6-d3', 'pm6-dh+', 'pm6-dh2', 'pm6-dh2x', 'pm6-d3h4', 'pm6-d3h4x', 'pm7', 'pm7-ts']

MULT_STRS = {1: 'SINGLET', 2: 'DOUBLET', 3: 'TRIPLET', 4: 'QUARTET', 5: 'QUINTET', 6: 'SEXTET', 7: 'SEPTET', 8: 'OCTET'}

base_cmd

Do only SCF AUX: Creates a checkpoint file NOREO: Dont reorient geometry

Type

1SCF

conf_key = 'mopac'

get_energy(atoms, coords)

Meant to be extended.

get_forces(atoms, coords)

Meant to be extended.

get_hessian(atoms, coords)

Meant to be extended.

parse_energy(path)

static parse_energy_from_aux(inp, *args, **kwargs)

parse_grad(path)

parse_hessian(path)

static parse_hessian_from_aux(inp, *args, **kwargs)

prepare_coords(atoms, coords, opt=False)

Get 3d coords in Angstrom.

Reshape internal 1d coords to 3d and convert to Angstrom.

Parameters

• atoms (iterable) -- Atom descriptors (element symbols).

• coords (np.array, 1d) -- 1D-array holding coordinates in Bohr.

Returns

coords -- 3D-array holding coordinates in Angstrom.

Return type

np.array, 3d

prepare_input(atoms, coords, calc_type, opt=False)

Meant to be extended.

read_aux(path)

6.5.2. Psi4

class pysisyphus.calculators.Psi4.Psi4(method, basis, to_set=None, pcm='iefpcm', solvent=None, write_fchk=False, **kwargs)

Bases: Calculator

conf_key = 'psi4'

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_fchk_str()

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)

Meant to be extended.

parse_energy(path)

parse_grad(path)

parse_hessian(path)

prepare_input(atoms, coords, calc_type)

Meant to be extended.

run_calculation(atoms, coords, **prepare_kwargs)

6.5.3. QCEngine

6.5.4. XTB 6.x

class pysisyphus.calculators.XTB.OptResult(opt_geom, opt_log)

Bases: tuple

opt_geom

Alias for field number 0

opt_log

Alias for field number 1

class pysisyphus.calculators.XTB.XTB(gbsa='', alpb='', gfn=2, acc=1.0, etemp=None, retry_etemp=None, topo=None, topo_update=None, quiet=False, **kwargs)

Bases: Calculator

__init__(gbsa='', alpb='', gfn=2, acc=1.0, etemp=None, retry_etemp=None, topo=None, topo_update=None, quiet=False, **kwargs)

XTB calculator.

Wrapper for running energy, gradient and Hessian calculations by XTB.

Parameters

• gbsa (str, optional) -- Solvent for GBSA calculation, by default no solvent model is used.

• alpb (str, optional) -- Solvent for ALPB calculation, by default no solvent model is used.

• gfn (int or str, must be (0, 1, 2, or "ff")) -- Hamiltonian for the XTB calculation (GFN0, GFN1, GFN2, or GFNFF).

• acc (float, optional) -- Accuracy control of the calculation, the lower the tighter several numerical thresholds are chosen.

• topo (str, optional) -- Path the a GFNFF-topolgy file. As setting up the topology may take some time for sizable systems, it may be desired to reuse the file.

• topo_update (int) -- Integer controlling the update interval of the GFNFF topology update. If supplied, the topolgy will be recreated every N-th calculation.

• mem (int) -- Mememory per core in MB.

• quiet (bool, optional) -- Suppress creation of log files.

static check_termination(inp, *args, **kwargs)

conf_key = 'xtb'

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_mdrestart_str(coords, velocities)

coords and velocities have to given in au!

get_pal_env()

get_retry_args()

keep(path)

Backup calculation results.

Parameters

path (Path) -- Temporary directory of the calculation.

Returns

kept_fns -- Dictonary holding the filenames that were backed up. The keys correspond to the type of file.

Return type

dict

parse_charges(fn=None)

parse_charges_from_json(fn=None)

parse_energy(path)

parse_gradient(path)

parse_hessian(path)

parse_md(path)

parse_opt(path, keep_log=False)

parse_topo(path)

prepare_add_args(xcontrol=None)

prepare_coords(atoms, coords)

Get 3d coords in Angstrom.

Reshape internal 1d coords to 3d and convert to Angstrom.

Parameters

• atoms (iterable) -- Atom descriptors (element symbols).

• coords (np.array, 1d) -- 1D-array holding coordinates in Bohr.

Returns

coords -- 3D-array holding coordinates in Angstrom.

Return type

np.array, 3d

prepare_input(atoms, coords, calc_type, point_charges=None)

Meant to be extended.

reattach(last_calc_cycle)

run_calculation(atoms, coords, **prepare_kwargs)

run_md(atoms, coords, t, dt, velocities=None, dump=1)

Expecting t and dt in fs, even though xtb wants t in ps!

run_opt(atoms, coords, keep=True, keep_log=False)

run_topo(atoms, coords)

write_mdrestart(path, mdrestart_str)

6.5.5. Dalton

class pysisyphus.calculators.Dalton.Dalton(basis, method='hf', **kwargs)

Bases: Calculator

compute(mol, prop)

conf_key = 'dalton'

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

prepare_input(atoms, coords)

Meant to be extended.

6.5.6. OpenBabel

class pysisyphus.calculators.OBabel.OBabel(ff='gaff', mol=None, **kwargs)

Bases: Calculator

conv_dict = {'kcal/mol': 627.5094740630558, 'kj/mol': 2625.4996394798254}

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

setup(atoms, coords)

6.6. Meta (wrapping) Calculators

6.6.1. ExternalPotential

6.6.2. Restraint

# General input structure for restraints
calc:
 type: ext
 # Multiple potentials could be specified here as a list
 potentials:
   # Primitive internal coordinate restraint
   - type: restraint
     # List of restraints; could also be multiple restraints. Each restraint is given as
     # list of 2 or 3 items.
     #
     # The first item always specifies an internal coordinate,
     # whereas the second argument is a force constant (in atomic units; actual units
     # depend on the coordinate). Optionally a reference value (third argument) can be
     # given. If omitted, the initial coordinate value is used as reference value.
     restraints: [[[BOND, 0, 1], 10, 3.0]]
     # The commented out input below would restrain the bond at its initial value.
     #restraints: [[[BOND, 0, 1], 10]]
     # Multiple restraints are specified as given below.
     #restraints: [[[BOND, 0, 1], 10], [[BEND, 0, 1, 2], 1.0]]
calc:
 type: [actual calculator that is wrapped by ExternalPotential]

6.6.3. HarmonicSphere

6.6.4. LogFermi

6.6.5. AFIR

class pysisyphus.calculators.AFIR.AFIR(calculator, fragment_indices, gamma, rho=1, p=6, ignore_hydrogen=True, complete_fragments=True, dump=True, h5_fn='afir.h5', h5_group_name='afir', **kwargs)

Bases: Calculator

__init__(calculator, fragment_indices, gamma, rho=1, p=6, ignore_hydrogen=True, complete_fragments=True, dump=True, h5_fn='afir.h5', h5_group_name='afir', **kwargs)

Artifical Force Induced Reaction calculator.

Currently, there are no automated drivers to run large-scale AFIR calculations with many different initial orientations and/or increasing collision energy parameter γ. Nontheless, selected AFIR calculations can be carried out by hand. After convergence, artificial potential & forces, as well as real energies and forces can be plotted with 'pysisplot --afir'. The highest energy point along the AFIR path can then be selected for a subsequent TS-optimization, e.g. via 'pysistrj --get [index] optimzation.trj'.

Future versions of pysisyphus may provide drivers for more automatted AFIR calculations.

Parameters

• calculator (Calculator) -- Actual QC calculator that provides energies and its derivatives, that are modified by the AFIR calculator, e.g., ORCA or Psi4.

• fragment_indices (List[List[int]]) -- List of lists of integers, specifying the separate fragments. If the indices in theses lists don't comprise all atoms in the molecule, the reamining indices will be added as a separate fragment. If a AFIR calculation is carried out with 2 fragments and 'complete_fragments' is True (see below) it is enough to specify only the indices of one fragment, e.g., for a system of 10 atoms 'fragment_indices=[[0,1,2,3]]' is enough. The second system will be set up automatically with indices [4,5,6,7,8,9].

• gamma (Union[_SupportsArray[dtype], _NestedSequence[_SupportsArray[dtype]], bool, int, float, complex, str, bytes, _NestedSequence[Union[bool, int, float, complex, str, bytes]]]) -- Collision energy parameter γ in au. For 2 fragments it can be a single integer, while for > 2 fragments a list of gammas must be given, specifying the pair-wise collision energy parameters. For 3 fragments 3 gammas must be given [γ_01, γ_02, γ_12], for 4 fragments 6 gammas would be required [γ_01, γ_02, γ_03, γ_12, γ_13, γ_23] and so on.

• rho (Union[_SupportsArray[dtype], _NestedSequence[_SupportsArray[dtype]], bool, int, float, complex, str, bytes, _NestedSequence[Union[bool, int, float, complex, str, bytes]]], default: 1) -- Direction of the artificial force, either 1 or -1. The same comments as for gamma apply. For 2 fragments a single integer is enough, for > 2 fragments a list of rhos must be given (see above). For rho=1 fragments are pushed together, for rho=-1 fragments are pulled apart.

• p (int, default: 6) -- Exponent p used in the calculation of the weight function ω. Defaults to 6 and probably does not have to be changed.

• ignore_hydrogen (bool, default: True) -- Whether hydrogens are ignored in the calculation of the artificial force.

• complete_fragments (bool, default: True) -- Whether an incomplete specification in 'fragment_indices' is automatically completed.

• dump (bool, default: True) -- Whether an HDF5 file is created.

• h5_fn (str, default: 'afir.h5') -- Filename of the HDF5 file used for dumping.

• h5_group_name (str, default: 'afir') -- HDF5 group name used for dumping.

• **kwargs -- Keyword arguments passed to the Calculator baseclass.

afir_fd_hessian_wrapper(coords3d, afir_grad_func)

property charge

dump_h5(atoms, coords, results)

get_energy(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)

Meant to be extended.

init_h5_group(atoms, max_cycles=None)

log_fragments()

property mult

set_atoms_and_funcs(atoms, coords)

Initially atoms was also an argument to the constructor of AFIR. I removed it so creation becomes easier. The first time a calculation is requested with a proper atom set everything is set up (cov. radii, afir function and corresponding gradient). Afterwards there is only a check if atoms != None and it is expected that all functions are properly set.

fragment_indices can also be incomplete w.r.t. to the number of atoms. If the sum of the specified fragment atoms is less than the number of atoms present then all remaining unspecified atoms will be gathered in one fragment.

write_fragment_geoms(atoms, coords)

class pysisyphus.calculators.AFIR.AFIRPath(atoms, cart_coords, energies, forces, opt_is_converged)

Bases: object

atoms: tuple

cart_coords: ndarray

energies: ndarray

forces: ndarray

opt_is_converged: bool

pysisyphus.calculators.AFIR.afir_closure(fragment_indices, cov_radii, gamma, rho=1, p=6, prefactor=1.0, logger=None)

rho=1 pushes fragments together, rho=-1 pulls fragments apart.

pysisyphus.calculators.AFIR.get_data_model(atoms, max_cycles)

6.6.6. ONIOM

class pysisyphus.calculators.ONIOMv2.LayerCalc(models, total_size, parent_layer_calc=None)

Bases: object

property charge

do_parent(with_parent)

static energy_from_results(model_energies, parent_energy=None)

get_energy(atoms, coords, with_parent=True)

get_forces(atoms, coords, with_parent=True)

get_hessian(atoms, coords, with_parent=True)

property mult

run_calculations(atoms, coords, method)

class pysisyphus.calculators.ONIOMv2.Link(ind, parent_ind, atom, g)

Bases: tuple

atom

Alias for field number 2

g

Alias for field number 3

ind

Alias for field number 0

parent_ind

Alias for field number 1

class pysisyphus.calculators.ONIOMv2.Model(name, calc_level, calc, parent_name, parent_calc_level, parent_calc, atom_inds, parent_atom_inds, use_link_atoms=True)

Bases: object

as_calculator(cap=False)

as_geom(all_atoms, all_coords)

capped_atoms_coords(all_atoms, all_coords)

create_bond_vec_getters(atoms)

create_links(atoms, coords, debug=False)

get_energy(atoms, coords, point_charges=None, parent_correction=True, cap=True)

get_forces(atoms, coords, point_charges=None, parent_correction=True, cap=True)

get_hessian(atoms, coords, point_charges=None, parent_correction=True, cap=True)

get_jacobian()

get_sparse_jacobian()

log(message='')

parse_charges()

class pysisyphus.calculators.ONIOMv2.ModelDummyCalc(model, cap=False)

Bases: object

get_energy(atoms, coords)

get_forces(atoms, coords)

class pysisyphus.calculators.ONIOMv2.ONIOM(calcs, models, geom, layers=None, embedding='', real_key='real', use_link_atoms=True, *args, **kwargs)

Bases: Calculator

__init__(calcs, models, geom, layers=None, embedding='', real_key='real', use_link_atoms=True, *args, **kwargs)

layer: list of models

len(layer) == 1: normal ONIOM, len(layer) >= 1: multicenter ONIOM.

model:

(sub)set of all atoms that resides in a certain layer and has a certain calculator.

atom_inds_in_layer(index, exclude_inner=False)

Returns list of atom indices in layer at index.

Atoms that also appear in inner layer can be excluded on request.

Parameters

• index (int) -- pasd

• exclude_inner (bool, default=False, optional) -- Whether to exclude atom indices that also appear in inner layers.

Returns

atom_indices -- List containing the atom indices in the selected layer.

Return type

list

calc_layer(atoms, coords, index, parent_correction=True)

property charge

embeddings = {'': '', 'electronic': 'Electronic embedding', 'electronic_rc': 'Electronic embedding with redistributed charges', 'electronic_rcd': 'Electronic embedding with redistributed charges and dipoles'}

get_energy(atoms, coords)

Meant to be extended.

get_forces(atoms, coords)

Meant to be extended.

get_hessian(atoms, coords)

Meant to be extended.

get_layer_calc(layer_ind)

property model_iter

property mult

run_calculation(atoms, coords)

run_calculations(atoms, coords, method)

pysisyphus.calculators.ONIOMv2.atom_inds_to_cart_inds(atom_inds)

pysisyphus.calculators.ONIOMv2.cap_fragment(atoms, coords, fragment, link_atom='H', g=0.709)

pysisyphus.calculators.ONIOMv2.get_embedding_charges(embedding, layer, parent_layer, coords3d)

pysisyphus.calculators.ONIOMv2.get_g_value(atom, parent_atom, link_atom)

6.6.7. Dimer

class pysisyphus.calculators.Dimer.Dimer(calculator, *args, N_raw=None, length=0.0189, rotation_max_cycles=15, rotation_method='fourier', rotation_thresh=0.0001, rotation_tol=1, rotation_max_element=0.001, rotation_interpolate=True, rotation_disable=False, rotation_disable_pos_curv=True, rotation_remove_trans=True, trans_force_f_perp=True, bonds=None, N_hessian=None, bias_rotation=False, bias_translation=False, bias_gaussian_dot=0.1, seed=None, write_orientations=True, forward_hessian=True, **kwargs)

Bases: Calculator

property C

Shortcut for the curvature.

property N

add_gaussian(atoms, center, N, height=0.1, std=0.0529, max_cycles=50, dot_ref=None)

property can_bias_f0

property can_bias_f1

property coords0

property coords1

curvature(f1, f2, N)

Curvature of the mode represented by the dimer.

direct_rotation(optimizer, prev_step)

do_dimer_rotations(rotation_thresh=None)

property energy0

property f0

property f1

property f1_bias

property f2

Never calculated explicitly, but estimated from f0 and f1.

fourier_rotation(optimizer, prev_step)

get_N_raw_from_hessian(h5_fn, root=0)

get_forces(atoms, coords)

Meant to be extended.

get_gaussian_energies(coords, sum_=True)

get_gaussian_forces(coords, sum_=True)

get_hessian(atoms, coords)

Meant to be extended.

remove_translation(displacement)

property rot_force

rotate_coords1(rad, theta)

Rotate dimer and produce new coords1.

set_N_raw(coords)

property should_bias_f0

May lead to calculation of f0 and/or f1 if present!

property should_bias_f1

May lead to calculation of f0 and/or f1 if present!

update_orientation(coords)

class pysisyphus.calculators.Dimer.Gaussian(height, center, std, N)

Bases: object

energy(R, height=None)

forces(R, height=None)

exception pysisyphus.calculators.Dimer.RotationConverged

Bases: Exception

6.7. Pure Python calculators

6.7.1. Sympy 2D Potentials

class pysisyphus.calculators.AnaPotBase.AnaPotBase(V_str, scale=1.0, xlim=(- 1, 1), ylim=(- 1, 1), levels=None, use_sympify=True, minima=None, saddles=None)

Bases: Calculator

anim_coords(coords, interval=50, show=False, title_func=None)

anim_opt(opt, energy_profile=False, colorbar=False, figsize=(8, 6), show=False)

get_energy(atoms, coords)

Meant to be extended.

get_forces(atoms, coords)

Meant to be extended.

classmethod get_geom(coords, atoms=('X',), V_str=None, calc_kwargs=None, geom_kwargs=None)

get_geoms_from_stored_coords(coords_list, i=None, calc_kwargs=None, geom_kwargs=None)

get_hessian(atoms, coords)

Meant to be extended.

get_minima(i=None, calc_kwargs=None, geom_kwargs=None)

get_path(num, minima_inds=None)

get_saddles(i=None, calc_kwargs=None, geom_kwargs=None)

plot(levels=None, show=False, **figkwargs)

plot3d(levels=None, show=False, zlim=None, vmin=None, vmax=None, resolution=100, rcount=50, ccount=50, nan_above=None, init_view=None, colorbar=False, **figkwargs)

plot_coords(xs, ys, enum=True, show=False, title=None)

plot_eigenvalue_structure(grid=50, levels=None, show=False)

plot_geoms(geoms, **kwargs)

plot_irc(irc, *args, **kwargs)

plot_opt(opt, *args, **kwargs)

statistics()

6.7.2. Lennard-Jones

class pysisyphus.calculators.LennardJones.LennardJones(sigma=1.8897261251, epsilon=1, rc=None)

Bases: Calculator

calculate(coords3d)

get_energy(atoms, coords)

Meant to be extended.

get_forces(atoms, coords)

Meant to be extended.

6.7.3. FakeASE

class pysisyphus.calculators.FakeASE.FakeASE(calc)

Bases: object

get_atoms_coords(atoms)

get_forces(atoms=None)

get_potential_energy(atoms=None)

6.7.4. TIP3P

class pysisyphus.calculators.TIP3P.TIP3P(rc=9.44863062728914)

Bases: Calculator

Transferable Intermolecular Potential 3 Point

aHOH = 104.52

calculate(coords3d)

charges

coulomb_energy = (multiple of elem. charge * multiple of elem. charge)

/ (distance in bohr) * 1 / (4 * pi * vacuum permittivity)

coulomb_prefactor converts everything to atomic units and it is ... drum roll ... 1. from scipy.constants import value as pcval self.coulomb_prefactor = (1 / (4 * np.pi) * pcval("elementary charge")**2

/ pcval("Hartree energy") / pcval("Bohr radius") / pcval("vacuum electric permittivity")

)

coulomb(coords3d)

epsilon = 0.0002423919586315716

get_energy(atoms, coords)

Meant to be extended.

get_forces(atoms, coords)

Meant to be extended.

qH = 0.417

qO = -0.834

rOH = 1.8088458464917874

sigma = 5.953790025507198