5. Coordinate Systems
Choice of coordinate system is often crucial, for the successful outcome of geometry optimizations. Pysisyphus supports different coordinate systems. By default, molecular optimizations are carried out in redundant internal coordinates (RIC).
5.1. Supported Coordinate Systems
5.1.1. Cartesian Coordinates
Strong coupling
Unambiguously defined, if translation and rotation (TR) are removed
Redundant set, if TR are not removed
type: cart
5.1.2. Redundant Internal Coordinates
Comprising bond stretches, (linear) bends, (improper) torsions and other primitives
Less coupling, compared to Cartesians
Easy estimation of diagonal model Hessians
Support constraints
Require sophisticated setup algorithm
Iterative internal-Cartesian backtransformation, which may fail
Usually highly redundant set
type: redund
5.1.3. Delocalized Internal Coordinates (DLC)
Complicated linear combinations of primitive internals
No coupling, coordinates are orthogonal to each other, at least at the geometrey they were defined at
Non redundant set
More efficient compared to RIC for bigger systems (if initial DLC generation is feasible)
Same comments apply, as for RICs
type: dlc
5.1.4. Translation & Rotation Internal Coordinates (TRIC)
Especially suited to optimize solvated/non-covalently bound systems
Translation and rotation coordinates are assigned to every fragment
Avoids error-prone assignment of interfragment coordinates
See 10.1063/1.4952956 for a full discussion
By default the B-Matrix is recalculated in every step of the internal-Cartesian backtransformation when TRIC is enabled
type: tric
5.2. Supported File Formats
All formats can be read, at least to a certain extent. Coordinate files are expected to be in Å, if not otherwise dictated by the respective format.
Suffix |
Write |
Comment |
---|---|---|
.xyz |
✓ |
Plain XYZ, single geometry. |
.trj |
✓ |
Plain XYZ, multiple geometries. |
.pdb |
✓ |
Protein Data Bank. |
.molden |
✗ |
Restricted to [Geometries] block. |
.zmat |
✗ |
Z-Matrix, see below for an example. |
.cjson |
✗ |
As saved by Avogadro. |
.crd |
✗ |
CHARMM card format |
.sdf |
✗ |
Structure-data file |
5.2.1. Z-Matrix example
C
N 1 1.35
H 1 1.0 2 105.0
H 2 1.4 1 105.0 3 150.0
H 2 1.4 1 110.0 3 -160.0
Indexing in Z-matrices is 1-based, so the first atom has index 1. Variable substitution is not supported in the current Z-matrix parser.
5.3. YAML Input
See below for an explanation of possible inputs in the geom section. Input related to internal coordinates is given mostly in the coord_kwargs subgroup, whereas atom-related input (isotops, freeze_atoms) is given one level above. See the example below:
geom:
type: cart # Coordinate system (cart/redund/dlc/tric)
fn: [input] # File name or inline input
union: False # Define same set of primitives at multiple geometries
isotopes: null # Specify different isotopes
freeze_atoms: null # Freeze Cartesians of certain atoms
coord_kwargs: # Keywords that are passed to the internal coordinate class
define_prims: null # Additionally define these primitives
constrain_prims: null # Primitive internals to be constrained
freeze_atoms_exclude: False # Whether to set up internal coordinates for frozen atoms
preopt:
geom: # geom block in preopt takes same keywords as above
... # no 'fn' key here!
endopt:
geom: # geom block in endopt takes same keywords as above
... # no 'fn' key here!
Employed coordinates and coordinate systems in preopt and endopt are similary controlled by a geom block. Same keywords are supported, as for the geom block, at the top level, except the fn key.
5.3.1. Inline input
Inline coordinates are supported for XYZ format, and expected in Å. Take care of proper indentation. The example below would yield RIC for the hydrogen molecule (1 bond).
geom:
type: redund
fn: |
2
H 0.0 0.0 0.0
H 0.0 0.0 0.7
5.3.2. Types of Primitive Coordinates
Pysisyphus implements many different (primitive) internal coordinates. Every coordinate is defined by its type and a set of atom indices, e.g., 2 indices for a bond, 3 indices for a bend and 4 indices for a dihedral.
Specification of a type is necessary, as there are many different kinds of bonds, bends and dihedrals/out-of-plane. One can't just assume, that a coordinate comprised of 3 atom indices is always a regular bend, as it may also be a linear bend or a translational coordinate (TRANSLATION_X, 14), describin the mean Cartesian X coordinate of 3 atoms.
Atom indices start at 0!
# Primitive types
BOND = 0
AUX_BOND = 1
HYDROGEN_BOND = 2
INTERFRAG_BOND = 3
AUX_INTERFRAG_BOND = 4
BEND = 5
LINEAR_BEND = 6
LINEAR_BEND_COMPLEMENT = 7
PROPER_DIHEDRAL = 8
IMPROPER_DIHEDRAL = 9
OUT_OF_PLANE = 10
LINEAR_DISPLACEMENT = 11
LINEAR_DISPLACEMENT_COMPLEMENT = 12
# TRANSLATION = 13 # Dummy coordinate
TRANSLATION_X = 14
TRANSLATION_Y = 15
TRANSLATION_Z = 16
# ROTATION = 17 # Dummy coordinate
ROTATION_A = 18
ROTATION_B = 19
ROTATION_C = 20
# CARTESIAN = 21 # Dummy coordinate
CARTESIAN_X = 22
CARTESIAN_Y = 23
CARTESIAN_Z = 24
BONDED_FRAGMENT = 25
DUMMY_TORSION = 26
DISTANCE_FUNCTION = 27
# atan2 based coordinates
BEND2 = 28
TORSION2 = 29
As some of these types are quite unwieldy, several shortcuts are supported, that can be used in place of the types above.
# Additional shortcuts
# Using Cartesians in the framework of internal coordinates is mainly
# useful if one wants to constrain certain atoms.
"X": [PT.CARTESIAN_X],
"Y": [PT.CARTESIAN_Y],
"Z": [PT.CARTESIAN_Z],
"XY": [PT.CARTESIAN_X, PT.CARTESIAN_Y],
"XZ": [PT.CARTESIAN_X, PT.CARTESIAN_Z],
"YZ": [PT.CARTESIAN_Y, PT.CARTESIAN_Z],
"XYZ": [PT.CARTESIAN_X, PT.CARTESIAN_Y, PT.CARTESIAN_Z],
"ATOM": [PT.CARTESIAN_X, PT.CARTESIAN_Y, PT.CARTESIAN_Z],
# Primitive aliases
"B": [PT.BOND],
"A": [PT.BEND],
"A2": [PT.BEND2],
"D": [PT.PROPER_DIHEDRAL],
"DIHEDRAL": [PT.PROPER_DIHEDRAL],
"TORSION": [PT.PROPER_DIHEDRAL],
"D2": [PT.PROPER_DIHEDRAL2],
"DIHEDRAL2": [PT.PROPER_DIHEDRAL2],
"TORSION2": [PT.PROPER_DIHEDRAL2],
# Translation & Rotation coordinates
"TRANSLATION": [PT.TRANSLATION_X, PT.TRANSLATION_Y, PT.TRANSLATION_Z],
"ROTATION": [PT.ROTATION_A, PT.ROTATION_B, PT.ROTATION_C],
# Miscellaneous
"DIST_FUNC": [PT.DISTANCE_FUNCTION]
5.3.3. Define Additional Primitives
Pysisyphus tries its best, to automatically come up with a reasonable set of internal coordinates, but sometimes the algorithm misses an important one. Especially at transition state guesses, where increased atom distances are common, bonds may be missed.
In such cases, additional coordinates can be requested explicitly. If additional coordinates are requested, a nested list is expected [[coord0], [coord1], ...].
# General structure (list of coordinate lists)
define_prims: [[PrimType or Shortcut], *[atom indices], ...]
# Examples
# Additional bond between atoms 4 and 7 (0-based indexing).
# All three lines below result in the same bond; the latter two use shortcuts.
define_prims: [[0, 4, 7]]
define_prims: [[B, 4, 7]]
define_prims: [[BOND, 4, 7]]
# Wrong specification (forgot outer list/brackets):
define_prims: [0, 4, 7]
# Also define an additional dihedral, beside the bond
define_prims: [[0, 4, 7], ["D", 0, 1, 2, 3]]
5.3.4. Freeze Atoms
All three Cartesian coordinates (X, Y, Z) of certain atoms can be frozen, so they always remain at their initial value. By setting freeze_atoms_exclude in coord_kwargs, frozen atoms can be excluded from the internal coordinate setup. By default frozen atoms are included.
geom:
type: [type]
fn: [fn]
freeze_atoms: [*atom indices]
coord_kwargs:
freeze_atoms_exclude: False
# Example; fully freeze Cartesians of first and second atom.
geom:
type: cart
fn: input.xyz
freeze_atoms: [0, 1]
coord_kwargs:
freeze_atoms: True
5.3.5. Constraints
Constraints beyond frozen atoms are currently only supported in conjunction with RIC (`type: redund`). It is not (yet) possible to modify the value of the specified coordinate via YAML input; the internal coordinate is constrained at its initial value. The same syntax as for define_prims is used. If the coordinate of the requested constraint is not already defined, it will be defined subsequently. There is no need to also add the constrained coordinate to define_prims.
# General structure (nested list of coordinates)
constrain_prims: [[[PrimType or Shortcut], *[atom indices]], ...]
# Examples
# Constrain Cartesian coordinate of atom 0.
# Both lines result in the same constraint.
constrain_prims: [[XYZ, 0]]
constrain_prims: [[ATOM, 0]]
# Constrain only Cartesian X and Y component of atom 0.
constrain_prims: [[XY, 0]]
# Constraint bond between atoms 4 and 7 (0-based indexing).
# All three lines below result in the same constraint; the latter two use shortcuts.
constrain_prims: [[0, 4, 7]]
constrain_prims: [[B, 4, 7]]
constrain_prims: [[BOND, 4, 7]]
Constraining the Cartesian coordinates (X, Y and Z) of one atom does not affect the final energy of an optimization. But constraining more than one atome does.
Harmonic restraints to selected primitive internals can be specified in the calc: section (see the Restraint documentation).
5.3.6. Isotopes
Different isotope masses can be requested. The system works similar to Gaussians system. A list of pairs is expected, where the first number specifies the atom and the second number is either an integer or a float. If it is an integer, the isotope mass closest to this integer is looked up in an internal database. Floats are used as is.
# General structure (nested list of coordinates)
isotopes: [[[atom index], [new mass, integer/float], ...]
# Modify the mass of atom with index 2 (hydrogen in this case)
# Both lines give identical results (deuterium).
# In the second line, the mass is given directly.
isotopes: [[2, 2]]
isotopes: [[2, 2.014101778]]
Different isotope masses affect calculated frequencies and IRCs. Atoms can be fixed in IRC calculations by specifying a very high mass. YAML does not recognize 1e9 as float, take care to add a dot (1.e9).
# Fix atom 0 in IRC calculation.
isotopes: [[0, 1.e9]]
5.4. Geometry & RedundantCoords
The central class in pysisyphus, handling coordinates and delegating calculations to (external QC) codes, is the Geometry class, similar to ASE's Atoms.
- class pysisyphus.Geometry.Geometry(atoms, coords, fragments=None, coord_type='cart', coord_kwargs=None, isotopes=None, freeze_atoms=None, remove_com=False, remove_centroid=False, comment='', name='')[source]
- __init__(atoms, coords, fragments=None, coord_type='cart', coord_kwargs=None, isotopes=None, freeze_atoms=None, remove_com=False, remove_centroid=False, comment='', name='')[source]
Object representing atoms in a coordinate system.
The Geometry represents atoms and their positions in coordinate system. By default cartesian coordinates are used, but internal coordinates are also possible.
- Parameters:
atoms (iterable) -- Iterable of length N, containing element symbols.
coords (1d iterable) -- 1d iterable of length 3N, containing the cartesian coordinates of N atoms.
fragments (dict, optional) -- Dict with different keys denoting different fragments. The values contain lists of atom indices.
coord_type ({"cart", "redund"}, optional) -- Type of coordinate system to use. Right now cartesian (cart) and redundand (redund) are supported.
coord_kwargs (dict, optional) -- Dictionary containing additional arguments that get passed to the constructor of the internal coordinate class.
isotopes (iterable of pairs, optional) -- Iterable of pairs consisting of 0-based atom index and either an integer or a float. If an integer is given the closest isotope mass will be selected. Given a float, this float will be directly used as mass.
freeze_atoms (iterable of integers) -- Specifies which atoms should remain fixed at their initial positions.
remove_com (bool, optional) -- Move center of mass to the origin.
remove_centroid (bool, optional) -- Move centroid to the origin.
comment (str, optional) -- Comment string.
name (str, optional) -- Verbose name of the geometry, e.g. methanal or water. Used for printing
- property all_energies
Return energies of all states that were calculated.
This will also set self.energy, which may NOT be the ground state, but the state correspondig to the 'root' attribute of the calculator.
- approximate_radius()[source]
Approximate molecule radius from the biggest atom distance along an axis.
- as_ascii_art()[source]
ASCII-art representation of the Geometry.
Using code from gpaw. Requires an ase installation.
- Return type:
str
- as_g98_list()[source]
Returns data for fake Gaussian98 standard orientation output.
- Returns:
g98_list -- List with one row per atom. Every row contains [center number, atomic number, atomic type (always 0 for now), X Y Z coordinates in Angstrom.
- Return type:
list
- as_xyz(comment='', atoms=None, cart_coords=None)[source]
Current geometry as a string in XYZ-format.
- Parameters:
comment (str, optional) -- Will be written in the second line (comment line) of the XYZ-string.
cart_coords (np.array, 1d, shape (3 * atoms.size, )) -- Cartesians for dumping instead of self._coords.
- Returns:
xyz_str -- Current geometry as string in XYZ-format.
- Return type:
str
- assert_compatibility(other)[source]
Assert that two Geometries can be substracted from each other.
- Parameters:
other (Geometry) -- Geometry for comparison.
- atom_indices()[source]
Dict with atom types as key and corresponding indices as values.
- Returns:
inds_dict -- Unique atom types as keys, corresponding indices as values.
- Return type:
dict
- property atom_types
- property atomic_numbers
- property bond_sets
- calc_relaxed_density(root, **prepare_kwargs)[source]
Calculate a relaxed excited state density via an ES gradient calculation.
The question is, if this method should set the wavefunction property at the current Geometry. On one hand, staying in pure python w/o numba the wavefunction sanity-check can become costly, even though it shouldn't be. On the other hand, setting the wavefunction would ensure consistency between the levels of theory used for density and wavefunction.
For now, calculating an ES density does not set a wavefunction on the Geometry, whereas requesting the relaxed density for the GS does.
TODO: add flag that allows setting the wavefunction (WF). Then, calculators should also include the WF in their results.
- property cart_coords
- property cart_forces
- property cart_gradient
- property cart_hessian
- property center_of_mass
Returns the center of mass.
- Returns:
R -- Center of mass.
- Return type:
np.array, shape(3, )
- center_of_mass_at(coords3d)[source]
Returns the center of mass at given coords3d.
- Parameters:
coords3d (np.array, shape(N, 3)) -- Cartesian coordiantes.
- Returns:
R -- Center of mass.
- Return type:
np.array, shape(3, )
- property centroid
Geometric center of the Geometry.
- Returns:
R -- Geometric center of the Geometry.
- Return type:
np.array, shape(3, )
- property comment
- coord_types = {'cart': None, 'cartesian': <class 'pysisyphus.intcoords.CartesianCoords.CartesianCoords'>, 'dlc': <class 'pysisyphus.intcoords.DLC.DLC'>, 'hdlc': <class 'pysisyphus.intcoords.DLC.HDLC'>, 'hredund': <class 'pysisyphus.intcoords.RedundantCoords.HybridRedundantCoords'>, 'mwcartesian': <class 'pysisyphus.intcoords.CartesianCoords.MWCartesianCoords'>, 'redund': <class 'pysisyphus.intcoords.RedundantCoords.RedundantCoords'>, 'tmtric': <class 'pysisyphus.intcoords.RedundantCoords.TMTRIC'>, 'tric': <class 'pysisyphus.intcoords.RedundantCoords.TRIC'>}
- property coords
1d vector of atomic coordinates.
- Returns:
coords -- 1d array holding the current coordinates.
- Return type:
np.array
- property coords3d
Coordinates in 3d.
- Returns:
coords3d -- Coordinates of the Geometry as 2D array.
- Return type:
np.array
- property coords_by_type
Coordinates in 3d by atom type and their corresponding indices.
- Returns:
cbt (dict) -- Dictionary with the unique atom types of the Geometry as keys. It's values are the 3d coordinates of the corresponding atom type.
inds (dict) -- Dictionary with the unique atom types of the Geometry as keys. It's values are the original indices of the 3d coordinates in the whole coords3d array.
- copy(coord_type=None, coord_kwargs=None)[source]
Returns a new Geometry object with same atoms and coordinates.
- Parameters:
coord_type (str) -- Desired coord_type, defaults to current coord_type.
coord_kwargs (dict, optional) -- Any desired coord_kwargs that will be passed to the RedundantCoords object.
- Returns:
geom -- New Geometry object with the same atoms and coordinates.
- Return type:
- property covalent_radii
- property energy
Energy of the current atomic configuration.
- Returns:
energy -- Energy of the current atomic configuration.
- Return type:
float
- property forces
Energy of the current atomic configuration.
- Returns:
force -- 1d array containing the forces acting on the atoms. Negative of the gradient.
- Return type:
np.array
- get_energy_and_forces_at(coords)[source]
Calculate forces and energies at the given coordinates.
The results are not saved in the Geometry object.
- get_normal_modes(cart_hessian=None, cart_gradient=None, proj_gradient=False, full=False)[source]
Normal mode wavenumbers, eigenvalues and Cartesian displacements Hessian.
- get_subgeom(indices, coord_type='cart', sort=False, cart_coords=None)[source]
Return a Geometry containing a subset of the current Geometry.
- Parameters:
indices (iterable of ints) -- Atomic indices that the define the subset of the current Geometry.
coord_type (str, ("cart", "redund"), optional) -- Coordinate system of the new Geometry.
cart_coords (default:
None
) -- Optional 1d array of Cartesian coordinates of shape (3*natoms, ).
- Returns:
sub_geom -- Subset of the current Geometry.
- Return type:
- property gradient
Negative of the force.
- Returns:
gradient -- 1d array containing the negative of the current forces.
- Return type:
np.array
- property has_energy
- property has_forces
- property hessian
Matrix of second derivatives of the energy in respect to atomic displacements.
- Returns:
hessian -- 2d array containing the second derivatives of the energy with respect to atomic/coordinate displacements depending on the type of coordiante system.
- Return type:
np.array
- property inertia_tensor
- property is_analytical_2d
- property is_linear
- jmol(atoms=None, cart_coords=None, stdin=None, jmol_cmd='jmol')[source]
Show geometry in jmol.
TODO: read jmol command from .pysisyphusrc ?!
- property layers
- property masses: ndarray
- property masses_rep
- property mm_inv
Inverted mass matrix.
Returns a diagonal matrix containing the inverted atomic masses.
- property mm_sqrt_inv
Inverted square root of the mass matrix.
- property moving_atoms
- property mw_coords
Mass-weighted coordinates.
- Returns:
mw_coords -- 1d array containing the mass-weighted cartesian coordiantes.
- Return type:
np.array
- property mw_gradient
Mass-weighted gradient.
- Returns:
mw_gradient -- Returns the mass-weighted gradient.
- Return type:
np.array
- property mw_hessian
Mass-weighted hessian.
- Returns:
mw_hessian -- 2d array containing the mass-weighted hessian M^(-1/2) H M^(-1/2).
- Return type:
np.array
- principal_axes_are_aligned()[source]
Check if the principal axes are aligned with the cartesian axes.
- Returns:
aligned -- Wether the principal axes are aligned or not.
- Return type:
bool
- set_coord(ind, coord)[source]
Set a coordinate by index.
- Parameters:
ind (int) -- Index in of the coordinate to set in the self.coords array.
coord (float) -- Coordinate value.
- set_results(results)[source]
Save the results from a dictionary.
- Parameters:
results (dict) -- The keys in this dict will be set as attributes in the current object, with the corresponding item as value.
- property sum_formula
- property td_1tdms
1-particle transition density matrices from TD-DFT/TDA.
Returns list of Xa, Ya, Xb and Yb in MO basis.
- property total_mass
- unweight_mw_hessian(mw_hessian)[source]
Unweight a mass-weighted hessian.
- Parameters:
mw_hessian (np.array) -- Mass-weighted hessian to be unweighted.
- Returns:
hessian -- 2d array containing the hessian.
- Return type:
np.array
- property vdw_radii
- property wavefunction
Basic infrastructure to work with internal coordinates is provided by RedundantCoords.
- class pysisyphus.intcoords.RedundantCoords.RedundantCoords(atoms, coords3d, masses=None, bond_factor=1.3, typed_prims=None, define_prims=None, constrain_prims=None, freeze_atoms=None, freeze_atoms_exclude=False, internals_with_frozen=False, define_for=None, bonds_only=False, check_bends=True, rebuild=True, bend_min_deg=15, dihed_max_deg=175, lb_min_deg=175, weighted=False, min_weight=0.3, svd_inv_thresh=0.000316, recalc_B=False, tric=False, hybrid=False, hbond_angles=False, rm_for_frag=None)[source]
- property B
Wilson B-Matrix
- property B_inv
Generalized inverse of the Wilson B-Matrix.
- property B_inv_prim
Generalized inverse of the primitive Wilson B-Matrix.
- property B_prim
Wilson B-Matrix
- property Bt_inv
Transposed generalized inverse of the Wilson B-Matrix.
- property Bt_inv_prim
Transposed generalized inverse of the primitive Wilson B-Matrix.
- property C
Diagonal matrix. Entries for constraints are set to one.
- property P
Projection matrix onto B. See [1] Eq. (4).
- backtransform_hessian(redund_hessian, int_gradient=None)[source]
Transform Hessian in internal coordinates to Cartesians.
- property bend_atom_indices
- property bend_indices
- property bond_atom_indices
- property bond_indices
- property bond_typed_prims
- property cartesian_indices
- property constrained_indices
- property coords
- property coords3d
- property dihedral_atom_indices
- property dihedral_indices
- property linear_bend_indices
- property outofplane_indices
- property prim_coords
- property prim_indices_set
- property prim_internals
- property primitives
- project_hessian(H, shift=1000)[source]
Expects a hessian in internal coordinates. See Eq. (11) in [1].
- property rotation_indices
- transform_hessian(cart_hessian, int_gradient=None)[source]
Transform Cartesian Hessian to internal coordinates.
- property translation_indices
- property typed_prims