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, keep_kind='all', check_mem=True, retry_calc=0, last_calc_cycle=None, clean_after=True, out_dir='qm_calcs', force_num_hess=False, num_hess_kwargs=None)[source]
__init__(calc_number=0, charge=0, mult=1, base_name='calculator', pal=1, mem=1000, keep_kind='all', check_mem=True, retry_calc=0, last_calc_cycle=None, clean_after=True, out_dir='qm_calcs', force_num_hess=False, num_hess_kwargs=None)[source]

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.

  • force_hess_kwargs (bool, default False) -- Force numerical Hessians.

  • num_hess_kwargs (dict) -- Keyword arguments for finite difference Hessian calculation.

apply_keep_kind()[source]
apply_set_plans(kept_fns, set_plans=None)[source]
build_set_plans(_set_plans=None)[source]
clean(path)[source]

Delete the temporary directory.

Parameters:

path (Path) -- Directory to delete.

conf_key = None
force_num_hessian()[source]

Always calculate numerical Hessians.

classmethod geom_from_fn(fn, **kwargs)[source]
get_cmd(key='cmd')[source]
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

get_num_hessian(atoms, coords, **prepare_kwargs)[source]
get_relaxed_density(atoms, coords, root, **prepare_kwargs)[source]

Meant to be extended.

get_restart_info()[source]

Return a dict containing chkfiles.

Returns:

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

Return type:

dict

get_stored_wavefunction(**kwargs)[source]
get_wavefunction(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

keep(path)[source]

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

load_wavefunction_from_file(fn, **kwargs)[source]
log(message='')[source]

Write a log message.

Wraps the logger variable.

Parameters:

message (str) -- Message to be logged.

make_fn(name, counter=None, return_str=False)[source]

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
popen(cmd, cwd=None)[source]
prepare(inp)[source]

Prepare a temporary directory and write input.

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

6. Paramters

inpstr

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

returns:
path: Path

Prepared directory.

prepare_coords(atoms, coords)[source]

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)[source]

Meant to be extended.

prepare_path(use_in_run=False)[source]

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: Path

Prepared directory.

prepare_pattern(raw_pat)[source]

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)[source]

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)[source]

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()[source]
print_out_fn(path)[source]

Print calculation output.

Prints the output of a calculator after a calculation.

Parameters:

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

restore_org_hessian()[source]

Restore original 'get_hessian' method, which may also fallback to numerical Hessians, if not implemented.

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)[source]

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)[source]

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)[source]

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)[source]

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

class pysisyphus.calculators.Calculator.HessKind(value)

An enumeration.

NUMERICAL = 2
ORG = 1
class pysisyphus.calculators.Calculator.KeepKind(value)

An enumeration.

ALL = 1
LATEST = 2
NONE = 3
class pysisyphus.calculators.Calculator.SetPlan(key, name=None, condition=<function SetPlan.<lambda>>, fail=None)[source]
condition()
fail: Optional[Callable] = None
key: str
name: Optional[str] = None

6.3. OverlapCalculator base class

class pysisyphus.calculators.OverlapCalculator.GroundStateContext(calc)[source]

Bases: object

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, root=None, nroots=None, track=False, ovlp_type='tden', double_mol=False, ovlp_with='previous', adapt_args=(0.5, 0.3, 0.6), cdds=None, orient='', dump_fn='overlap_data.h5', h5_dump=False, conf_thresh=0.001, mos_ref='cur', mos_renorm=True, min_cost=False, **kwargs)[source]

Bases: Calculator

H5_MAP = {'Ca': 'Ca_list', 'Cb': 'Cb_list', 'Xa': 'Xa_list', 'Xb': 'Xb_list', 'Ya': 'Ya_list', 'Yb': 'Yb_list', 'all_energies': 'all_energies_list', 'coords': 'coords_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', 'top': 'transition orbital pair overlap', 'wf': 'wavefunction overlap'}
VALID_CDDS = (None, 'calc', 'render')
VALID_KEYS = ['wf', 'tden', 'nto', 'nto_org', 'top']
VALID_XY = ('X', 'X+Y', 'X-Y')
calc_cdd_cube(root, cycle=-1)[source]
conf_thresh

self.dyn_roots = int(dyn_roots) if self.dyn_roots != 0:

self.dyn_roots = 0 self.log("dyn_roots = 0 is hardcoded right now")

property data_model
dump_overlap_data()[source]
static from_overlap_data(h5_fn, set_wfow=False)[source]
get_ci_coeffs_for(ind)[source]
get_h5_group()[source]
get_indices(indices=None)[source]

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(C, S_AO)[source]
get_orbital_matrices(indices=None, S_AO=None)[source]

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(C_ref, C_cur)[source]
static get_sao_from_mo_coeffs(C)[source]

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

get_tden_overlaps(indices=None, S_AO=None)[source]
get_top_differences(indices=None, S_AO=None)[source]

Transition orbital projection.

prepare_overlap_data(path)[source]

This method has to implement the calculator specific parsing of MO-coefficients, CI-coefficients and energies. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

render_cdd_cube()[source]
static renorm_mos(C, S_AO)[source]
store_overlap_data(atoms, coords, path=None, overlap_data=None)[source]
property stored_calculations
track_root(ovlp_type=None)[source]

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

pysisyphus.calculators.OverlapCalculator.get_data_model(exc_state_num, occ_a, virt_a, occ_b, virt_b, ovlp_type, atoms, max_cycles)[source]

6.4. Calculators with Excited state capabilities

6.4.1. Gaussian09

class pysisyphus.calculators.Gaussian09.Gaussian09(*args, **kwargs)[source]
conf_key = 'gaussian09'

6.4.2. Gaussian16

class pysisyphus.calculators.Gaussian16.Gaussian16(route, gbs='', gen='', keep_chk=False, wavefunction_dump=True, stable='', fchk=None, iop9_40=3, **kwargs)[source]

Bases: OverlapCalculator

conf_key = 'gaussian16'
get_all_energies(atoms, coords, **prepare_kwargs)[source]
get_chkfiles()[source]
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

get_relaxed_density(atoms, coords, root, **prepare_kwargs)[source]

Meant to be extended.

get_stored_wavefunction(**kwargs)[source]
get_wavefunction(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

make_exc_str()[source]
make_fchk(path)[source]
make_gbs_str()[source]
parse_635r_dump(dump_path, roots, nmos)[source]
parse_all_energies(fchk=None)[source]
parse_charges(path=None)[source]
parse_double_mol(path, out_fn=None)[source]
parse_energy(path)[source]
static parse_fchk(fchk_path, keys)[source]
parse_force(path)[source]
parse_hessian(path)[source]
parse_keyword(text)[source]
parse_log(*args, **kwargs)
parse_stable(path)[source]
parse_tddft(path)[source]
prepare_input(atoms, coords, calc_type, did_stable=False, point_charges=None)[source]

Meant to be extended.

prepare_overlap_data(path)[source]

This method has to implement the calculator specific parsing of MO-coefficients, CI-coefficients and energies. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

reuse_data(path)[source]
run_after(path)[source]

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)[source]
run_double_mol_calculation(atoms, coords1, coords2)[source]
run_rwfdump(path, rwf_index, chk_path=None)[source]
run_stable(atoms, coords, **prepare_kwargs)[source]
set_chkfiles(chkfiles)[source]
store_and_track(results, func, atoms, coords, **prepare_kwargs)[source]
pysisyphus.calculators.Gaussian16.get_nmos(text)[source]
pysisyphus.calculators.Gaussian16.noparse(path)[source]
pysisyphus.calculators.Gaussian16.to_ind_and_spin(lbl)[source]

6.4.3. OpenMolcas

Pysisyphus currently supports energy and gradient calculations utilizing the &rasscf and/or the &mcpdft sections. Neither analytical nor numerical Hessians are yet implemented for the OpenMolcas-calculator.

Two keywords are always required: inporb and basis, with the former pointing to a .RasOrb file and the latter containing the selected atomic orbital basis, e.g., ano-rcc-vdzp. Additional input for the &gateway, &rasscf and &mcpdft sections can be given via the respective keyword(s).

Due to restrictions of the current design, simple keywords that don't take further arguments as cmsi in &rasscf or grad and mspdft in &mcpdft still must be given with a trailing colon. See below for an example.

Listing 6.1 Optimization using compressed-multi-state PDFT.
geom:
 type: redund
 fn: |
  4
  
  C	-1.0398336639	 0.0	    0.0
  S	 0.6002429216	 0.0	    0.0
  H -1.6321592382	-0.94139864	0.0
  H	-1.6321592382	 0.94139864	0.0
calc:
 type: openmolcas
 basis: cc-pvdz
 charge: 0
 mult: 1
 inporb: /home/johannes/tmp/359_cmspdft/359_cmspdft.RasOrb
 rasscf:
  ciroot: 3 3 1
  mdrlxroot: 2
  cmsi:
 mcpdft:
  ksdft: t:pbe
  grad:
  mspdft:
opt:
 thresh: gau
class pysisyphus.calculators.OpenMolcas.OpenMolcas(basis, inporb, rasscf=None, gateway=None, mcpdft=None, rassi=None, nprocs=1, track=True, omp_var='OMP_NUM_THREADS', **kwargs)[source]

Bases: Calculator

build_gateway_str()[source]
build_mcpdft_str()[source]
build_rasscf_str()[source]
build_rassi_str()[source]
build_str_from_dict(dct)[source]
conf_key = 'openmolcas'
get_energy(atoms, coords)[source]

Meant to be extended.

get_forces(atoms, coords)[source]

Meant to be extended.

get_pal_env()[source]
get_root()[source]
parse_energies(text)[source]
parse_energy(path)[source]
parse_gradient(path)[source]
parse_rassi_track(path)[source]
prepare_coords(atoms, coords)[source]

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)[source]

Meant to be extended.

reattach(last_calc_cycle)[source]
run_calculation(atoms, coords, calc_type='energy')[source]

6.4.4. ORCA 4.2.1 / 5.0.1

class pysisyphus.calculators.ORCA.ORCA(keywords, blocks='', gbw=None, do_stable=False, numfreq=False, json_dump=None, wavefunction_dump=True, **kwargs)[source]

Bases: OverlapCalculator

__init__(keywords, blocks='', gbw=None, do_stable=False, numfreq=False, json_dump=None, wavefunction_dump=True, **kwargs)[source]

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 (bool, optional) -- Use numerical frequencies instead of analytical ones.

  • json_dump (bool, optional, deprecated) -- Use 'wavefunction_dump' instead.

  • wavefunction_dump (bool, optional) -- Whether to dump the wavefunction to BSON via orca_2json. The BSON can become very large in calculations comprising many basis functions.

check_termination(*args, **kwargs)
clean_tmp(path)[source]
conf_key = 'orca'
get_all_energies(atoms, coords, **prepare_kwargs)[source]
get_block_str()[source]
get_chkfiles()[source]
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

get_moinp_str(gbw)[source]
get_relaxed_density(atoms, coords, root, **prepare_kwargs)[source]

Meant to be extended.

get_stable_wavefunction(atoms, coords)[source]
get_stored_wavefunction(**kwargs)[source]
get_wavefunction(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

parse_all_energies(text=None, triplets=None)[source]
parse_all_energies_from_path(path)[source]
static parse_atoms_coords(inp, *args, **kwargs)
static parse_cis(cis)[source]

Simple wrapper of external function.

Currently, only returns Xα and Yα.

parse_energy(path)[source]
parse_engrad(path)[source]
static parse_engrad_info(inp, *args, **kwargs)
static parse_gbw(gbw_fn)[source]
static parse_hess_file(inp, *args, **kwargs)
parse_hessian(path)[source]
parse_mo_numbers(out_fn)[source]
parse_stable(path)[source]
prepare_input(atoms, coords, calc_type, point_charges=None, do_stable=False)[source]

Meant to be extended.

prepare_overlap_data(path, triplets=None)[source]

This method has to implement the calculator specific parsing of MO-coefficients, CI-coefficients and energies. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

reattach(last_calc_cycle)[source]
run_after(path)[source]

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)[source]

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

set_chkfiles(chkfiles)[source]
set_mo_coeffs(mo_coeffs=None, gbw=None)[source]
static set_mo_coeffs_in_gbw(in_gbw_fn, out_gbw_fn, mo_coeffs)[source]

See self.parse_gbw.

store_and_track(results, func, atoms, coords, **prepare_kwargs)[source]
class pysisyphus.calculators.ORCA.ORCAGroundStateContext(calc)[source]

Bases: GroundStateContext

pysisyphus.calculators.ORCA.geom_from_orca_hess(fn)[source]
pysisyphus.calculators.ORCA.get_exc_ens_fosc(wf_fn, cis_fn, log_fn)[source]
pysisyphus.calculators.ORCA.get_name(text)[source]

Return string that comes before first character & offset.

pysisyphus.calculators.ORCA.make_sym_mat(table_block)[source]
pysisyphus.calculators.ORCA.parse_orca_cis(cis_fn, restricted_same_ab=False, triplets_only=True)[source]
Read binary CI vector file from ORCA.

Loosly based on TheoDORE 1.7.1, Authors: S. Mai, F. Plasser https://sourceforge.net/p/theodore-qc

With restricted_same_ab the alpha part will be copied over to the beta-part in restricted calculations. Otherwise the beta-part (Xb, Yb) will just be zeros in restricted calculations.

pysisyphus.calculators.ORCA.parse_orca_gbw(gbw_fn)[source]

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

pysisyphus.calculators.ORCA.parse_orca_gbw_new(gbw_fn)[source]

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

Return type:

MOCoeffs

pysisyphus.calculators.ORCA.save_orca_pc_file(point_charges, pc_fn, hardness=None)[source]
pysisyphus.calculators.ORCA.update_gbw(gbw_in, gbw_out, alpha_mo_coeffs=None, beta_mo_coeffs=None, alpha_energies=None, alpha_occs=None, beta_energies=None, beta_occs=None)[source]

MOs are expected to be in columns.

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.ExSpectrumRoot(root, sym, exc_energy, osc_vel, osc_len)[source]

Bases: object

exc_energy: float
osc_len: float
osc_vel: float
root: int
sym: str
class pysisyphus.calculators.Turbomole.Turbomole(control_path=None, numfreq=False, simple_input=None, double_mol_path=None, cosmo_kwargs=None, wavefunction_dump=True, **kwargs)[source]

Bases: OverlapCalculator

append_control(to_append, log_msg='', **kwargs)[source]
conf_key = 'turbomole'
get_all_energies(atoms, coords, **prepare_kwargs)[source]
get_chkfiles()[source]
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, cmd=None, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

get_pal_env()[source]
get_relaxed_density(atoms, coords, root, **prepare_kwargs)[source]

Meant to be extended.

get_ricc2_root(text)[source]
get_stored_wavefunction(**kwargs)[source]
get_wavefunction(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

static make_molden(path)[source]
parse_all_energies()[source]
parse_cc2_vectors(ccre)[source]
parse_ci_coeffs()[source]
parse_double_mol(path)[source]

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

parse_energy(path)[source]
parse_force(path)[source]
parse_gs_energy()[source]

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)[source]
parse_mos()[source]
static parse_td_vectors(text)[source]

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.

static parse_tddft_tden(inp, *args, **kwargs)
prepare_input(atoms, coords, calc_type, point_charges=None)[source]

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)[source]

This method has to implement the calculator specific parsing of MO-coefficients, CI-coefficients and energies. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

prepare_point_charges(point_charges)[source]

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

prepare_td(text)[source]
run_after(path)[source]

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)[source]
run_double_mol_calculation(atoms, coords1, coords2)[source]
set_chkfiles(chkfiles)[source]
set_occ_and_mo_nums(text)[source]
store_and_track(results, func, atoms, coords, **prepare_kwargs)[source]
sub_control(pattern, repl, log_msg='', **kwargs)[source]
class pysisyphus.calculators.Turbomole.TurbomoleGroundStateContext(calc)[source]

Bases: GroundStateContext

pysisyphus.calculators.Turbomole.control_from_simple_input(simple_inp, charge, mult, cosmo_kwargs=None)[source]

Create control file from 'simple input'.

See examples/opt/26_turbomole_simple_input/ for an example.

pysisyphus.calculators.Turbomole.get_cosmo_data_groups(atoms, epsilon, rsolv=None, refind=None, dcosmo_rs=None)[source]
pysisyphus.calculators.Turbomole.get_density_matrices_for_root(log_fn, vec_fn, root, rlx_vec_fn=None, Ca=None, Cb=None)[source]
pysisyphus.calculators.Turbomole.index_strs_for_atoms(atoms)[source]
pysisyphus.calculators.Turbomole.parse_frozen_nmos(text)[source]

Determine number of occ. & and virt. orbitals used in ES calculations.

Return type:

tuple[list[tuple[int, int], tuple[int, int]], bool]

pysisyphus.calculators.Turbomole.render_data_groups(raw_data_groups)[source]

6.4.7. DFTB+ 20.x

class pysisyphus.calculators.DFTBp.DFTBp(parameter, *args, slakos=None, **kwargs)[source]

Bases: OverlapCalculator

conf_key = 'dftbp'
get_all_energies(atoms, coords, **prepare_kwargs)[source]
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

static get_excited_state_str(track, root, nroots, forces=False)[source]
get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

static get_gen_str(atoms, coords)[source]
hubbard_derivs = {'3ob': {'Br': -0.0573, 'C': -0.1492, 'Ca': -0.034, 'Cl': -0.0697, 'F': -0.1623, 'H': -0.1857, 'I': -0.0433, 'K': -0.0339, 'Mg': -0.02, 'N': -0.1535, 'Na': -0.0454, 'O': -0.1575, 'P': -0.14, 'S': -0.11, 'Zn': -0.03}}
max_ang_moms = {'3ob': {'Br': 'd', 'C': 'p', 'Ca': 'p', 'Cl': 'd', 'F': 'p', 'H': 's', 'I': 'd', 'K': 'p', 'Mg': 'p', 'N': 'p', 'Na': 'p', 'O': 'p', 'P': 'd', 'S': 'd', 'Zn': 'd'}, 'mio-ext': {'C': 'p', 'H': 's', 'N': 'p', 'O': 'p'}}
parse_all_energies(out_fn=None, exc_dat=None)[source]
parse_energy(path)[source]
static parse_exc_dat(text)[source]
parse_forces(path)[source]
parse_total_energy(text)[source]
prepare_input(atoms, coords, calc_type)[source]

Meant to be extended.

prepare_overlap_data(path)[source]

This method has to implement the calculator specific parsing of MO-coefficients, CI-coefficients and energies. Should return a filename pointing to TURBOMOLE like mos, a MO coefficient array and a CI coefficient array.

run_calculation(atoms, coords, **prepare_kwargs)[source]
store_and_track(results, func, atoms, coords, **prepare_kwargs)[source]
pysisyphus.calculators.DFTBp.parse_mo(eigvec)[source]
pysisyphus.calculators.DFTBp.parse_xplusy(text)[source]

6.5. Calculators with Ground state capabilities

6.5.1. MOPAC 2016

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

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, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

parse_energy(path)[source]
static parse_energy_from_aux(inp, *args, **kwargs)
parse_grad(path)[source]
parse_hessian(path)[source]
static parse_hessian_from_aux(inp, *args, **kwargs)
prepare_coords(atoms, coords, opt=False)[source]

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)[source]

Meant to be extended.

read_aux(path)[source]
run_calculation(atoms, coords, **prepare_kwargs)[source]

6.5.2. Psi4

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

Bases: Calculator

conf_key = 'psi4'
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_fchk_str()[source]
get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

parse_energy(path)[source]
parse_grad(path)[source]
parse_hessian(path)[source]
prepare_input(atoms, coords, calc_type)[source]

Meant to be extended.

run_calculation(atoms, coords, **prepare_kwargs)[source]

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, restart=False, topo=None, topo_update=None, quiet=False, wavefunction_dump=False, **kwargs)[source]

Bases: Calculator

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

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.

  • wavefunction_dump (bool) -- Whether to dump a molden file.

static check_termination(inp, *args, **kwargs)
conf_key = 'xtb'
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

get_mdrestart_str(coords, velocities)[source]

coords and velocities have to given in au!

get_pal_env()[source]
get_retry_args()[source]
get_stored_wavefunction(**kwargs)[source]
parse_charges(fn=None)[source]
parse_charges_from_json(fn=None)[source]
parse_energy(path)[source]
parse_gradient(path)[source]
parse_hessian(path)[source]
parse_md(path)[source]
parse_opt(path, keep_log=False)[source]
parse_topo(path)[source]
prepare_add_args(xcontrol=None)[source]
prepare_coords(atoms, coords)[source]

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)[source]

Meant to be extended.

reattach(last_calc_cycle)[source]
run_calculation(atoms, coords, **prepare_kwargs)[source]
run_md(atoms, coords, t, dt, velocities=None, dump=1)[source]

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

run_opt(atoms, coords, keep=True, keep_log=False)[source]
run_topo(atoms, coords)[source]
write_mdrestart(path, mdrestart_str)[source]

6.5.5. Dalton

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

Bases: Calculator

compute(mol, prop)[source]
conf_key = 'dalton'
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

prepare_input(atoms, coords)[source]

Meant to be extended.

6.5.6. OpenBabel

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

Bases: Calculator

conv_dict = {'kcal/mol': 627.5094740630558, 'kj/mol': 2625.4996394798254}
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

setup(atoms, coords)[source]

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. RMSD

6.6.6. DFT-D3

Method to add DFT-D3 dispersion corrections as an external potential via the program developed by the Grimme group <https://www.chemie.uni-bonn.de/grimme/de/software/dft-d3/get_dft-d3>.

This is for use with calculators that do not natively provide D3 corrections (e.g. OpenMolcas). Usage mirrors that of other external potentials, with an example given below.

# General input structure for restraints
calc:
 type: ext
 # Multiple potentials could be specified here as a list
 potentials:
   # Add atom-pairwise D3 dispersion correction as a differentiable, external potential
   - type: d3
     # Functional is specified in TURBOMOLE format, all lower case.
     functional: pbe
     # Optional Becke-Johnson damping, default false, recommended true
     bjdamping: true

calc:
 type: [actual calculator that is wrapped by ExternalPotential]
class pysisyphus.calculators.DFTD3.DFTD3(geom, functional, bjdamping=False, **kwargs)[source]
calc(coords3d, gradient=False)[source]
conf_key = 'dftd3'
parse_energy(path)[source]
parse_gradient(path)[source]

6.6.7. AFIR

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

Bases: Calculator

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

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[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], 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[Any]], _NestedSequence[_SupportsArray[dtype[Any]]], 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: False) -- Whether hydrogens are ignored in the calculation of the artificial force. All weights between atom pairs containing hydrogen will be set to 0.

  • zero_hydrogen (bool, default: True) -- Whether to use 0.0 as covalent radius for hydrogen in the weight function. Compared to 'ignore_hydrogen', which results in zero weights for all atom pairs involving hydrogen, 'zero_hydrogen' may be non-zero, depending on the covalent radius of the second atom in the pair.

  • 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)[source]
property charge
dump_h5(atoms, coords, results)[source]
get_energy(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_forces(atoms, coords, **prepare_kwargs)[source]

Meant to be extended.

get_hessian(atoms, coords, **prepare_kwargs)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

init_h5_group(atoms, max_cycles=None)[source]
log_fragments()[source]
property mult
set_atoms_and_funcs(atoms, coords)[source]

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)[source]
exception pysisyphus.calculators.AFIR.CovRadiiSumZero[source]

Bases: Exception

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

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

pysisyphus.calculators.AFIR.get_data_model(atoms, max_cycles)[source]

6.6.8. ONIOM

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

Bases: object

property charge
do_parent(with_parent)[source]
static energy_from_results(model_energies, parent_energy=None)[source]
get_energy(atoms, coords, with_parent=True)[source]
get_forces(atoms, coords, with_parent=True)[source]
get_hessian(atoms, coords, with_parent=True)[source]
property mult
run_calculations(atoms, coords, method)[source]

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)[source]

Bases: object

as_calculator(cap=False)[source]
as_geom(all_atoms, all_coords)[source]
capped_atoms_coords(all_atoms, all_coords)[source]
create_bond_vec_getters(atoms)[source]
get_energy(atoms, coords, point_charges=None, parent_correction=True, cap=True)[source]
get_forces(atoms, coords, point_charges=None, parent_correction=True, cap=True)[source]
get_hessian(atoms, coords, point_charges=None, parent_correction=True, cap=True)[source]
get_jacobian()[source]
get_sparse_jacobian()[source]
log(message='')[source]
parse_charges()[source]
class pysisyphus.calculators.ONIOMv2.ModelDummyCalc(model, cap=False)[source]

Bases: object

get_energy(atoms, coords)[source]
get_forces(atoms, coords)[source]
class pysisyphus.calculators.ONIOMv2.ONIOM(calcs, models, geom, layers=None, embedding='', real_key='real', use_link_atoms=True, *args, **kwargs)[source]

Bases: Calculator

__init__(calcs, models, geom, layers=None, embedding='', real_key='real', use_link_atoms=True, *args, **kwargs)[source]
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)[source]

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)[source]
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)[source]

Meant to be extended.

get_forces(atoms, coords)[source]

Meant to be extended.

get_hessian(atoms, coords)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

get_layer_calc(layer_ind)[source]
property model_iter
property mult
run_calculation(atoms, coords)[source]
run_calculations(atoms, coords, method)[source]
pysisyphus.calculators.ONIOMv2.atom_inds_to_cart_inds(atom_inds)[source]
pysisyphus.calculators.ONIOMv2.cap_fragment(atoms, coords, fragment, link_atom='H', g=None)[source]
pysisyphus.calculators.ONIOMv2.get_embedding_charges(embedding, layer, parent_layer, coords3d)[source]
pysisyphus.calculators.ONIOMv2.get_g_value(atom, parent_atom, link_atom)[source]

6.6.9. 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)[source]

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)[source]
property can_bias_f0
property can_bias_f1
property coords0
property coords1
curvature(f1, f2, N)[source]

Curvature of the mode represented by the dimer.

direct_rotation(optimizer, prev_step)[source]
do_dimer_rotations(rotation_thresh=None)[source]
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)[source]
get_N_raw_from_hessian(h5_fn, root=0)[source]
get_forces(atoms, coords)[source]

Meant to be extended.

get_gaussian_energies(coords, sum_=True)[source]
get_gaussian_forces(coords, sum_=True)[source]
get_hessian(atoms, coords)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

remove_translation(displacement)[source]
property rot_force
rotate_coords1(rad, theta)[source]

Rotate dimer and produce new coords1.

set_N_raw(coords)[source]
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)[source]
class pysisyphus.calculators.Dimer.Gaussian(height, center, std, N)[source]

Bases: object

energy(R, height=None)[source]
forces(R, height=None)[source]
exception pysisyphus.calculators.Dimer.RotationConverged[source]

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, **kwargs)[source]

Bases: Calculator

anim_coords(coords, interval=50, show=False, title_func=None)[source]
anim_opt(opt, energy_profile=False, colorbar=False, figsize=(8, 6), show=False)[source]
get_energy(atoms, coords)[source]

Meant to be extended.

get_forces(atoms, coords)[source]

Meant to be extended.

classmethod get_geom(coords, atoms=('X',), V_str=None, calc_kwargs=None, geom_kwargs=None)[source]
get_geoms_from_stored_coords(coords_list, i=None, calc_kwargs=None, geom_kwargs=None)[source]
get_hessian(atoms, coords)[source]

Get Hessian matrix. Fall back to numerical Hessian, if not overriden.

Preferrably, this method should provide an analytical Hessian.

get_minima(i=None, calc_kwargs=None, geom_kwargs=None)[source]
get_path(num, minima_inds=None)[source]
get_saddles(i=None, calc_kwargs=None, geom_kwargs=None)[source]
plot(levels=None, show=False, **figkwargs)[source]
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)[source]
plot_coords(xs, ys, enum=True, show=False, title=None)[source]
plot_eigenvalue_structure(grid=50, levels=None, show=False)[source]
plot_geoms(geoms, **kwargs)[source]
plot_irc(irc, *args, **kwargs)[source]
plot_opt(opt, *args, **kwargs)[source]
statistics()[source]

6.7.2. Lennard-Jones

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

Bases: Calculator

calculate(coords3d)[source]
get_energy(atoms, coords)[source]

Meant to be extended.

get_forces(atoms, coords)[source]

Meant to be extended.

6.7.3. FakeASE

class pysisyphus.calculators.FakeASE.FakeASE(calc)[source]

Bases: object

get_atoms_coords(atoms)[source]
get_forces(atoms=None)[source]
get_potential_energy(atoms=None)[source]

6.7.4. TIP3P

class pysisyphus.calculators.TIP3P.TIP3P(rc=9.44863062728914)[source]

Bases: Calculator

Transferable Intermolecular Potential 3 Point

aHOH = 104.52
calculate(coords3d)[source]
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)[source]
epsilon = 0.0002423919586315716
get_energy(atoms, coords)[source]

Meant to be extended.

get_forces(atoms, coords)[source]

Meant to be extended.

qH = 0.417
qO = -0.834
rOH = 1.8088458464917874
sigma = 5.953790025507198