Source code for pysisyphus.drivers.birkholz

# [1]
#     Birkholz, Schlegel, 2015

import numpy as np

from pysisyphus.Geometry import Geometry
from pysisyphus.drivers import relaxed_scan
from pysisyphus.helpers_pure import highlight_text
from pysisyphus.intcoords.setup import get_bond_mat
from pysisyphus.intcoords.PrimTypes import normalize_prim_inputs

[docs] def bond_order(r, r0, b=2): """Bond order for given bond length and reference length. Eq. (3) in [1].""" return max(0, (b * (r0 / r) - 1) / (b - 1))
[docs] def bond_orders(coords3d, bond_indices, r0s): """List of bond orders.""" bos = list() for r0, (a, b) in zip(r0s, bond_indices): r = np.linalg.norm(coords3d[a] - coords3d[b]) bo = bond_order(r, r0) bos.append(bo) return np.array(bos)
[docs] def get_r0s(geom, bond_indices): """Reference bond lengths as sum of covalent radii.""" crs = geom.covalent_radii a, b = bond_indices.T return crs[a] + crs[b]
[docs] def bond_orders_for_geom(geom, bond_indices): """Wrapper for bond_orders for simple use with Geometry.""" r0s = get_r0s(geom, bond_indices) return bond_orders(geom.coords3d, bond_indices, r0s)
[docs] def length_for_bond_order(bo, r0, b=2): """Return bond length for given bond order and reference length. Eq. (3) in [1].""" return b / ((b - 1) * bo + 1) * r0
[docs] def birkholz_interpolation(geoms, calc_getter, recreate=True): assert len(geoms) >= 2 start = geoms[0] end = geoms[-1] geoms = (start, end) atoms = start.atoms print(f"Coordinates of 'start' geometry\n{start.as_xyz()})") print(f"Coordinates of 'end' geometry\n{end.as_xyz()})\n") print(highlight_text("Bonding analysis")) # Bond matrices bm_start, bm_end = [get_bond_mat(geom) for geom in geoms] # Determine formed/broken bonds. The lower triangular part of the bond # matrices has to be zeroed, to avoid double counting. Changes in bonding # are determined with logical XOR. # # Start End XOR Comment # ----- --- --- ------- # True True False Bond is present at both geometries. # False False False Bond is absent at both geometries. # True False True Bond is broken when going from start to end. # False True True Bond is formed when going from start to end. bm_diff = np.triu(bm_start) ^ np.triu(bm_end) bond_indices = np.stack(np.nonzero(bm_diff), axis=1) bonding_changes = bm_diff.sum() print(f"{bonding_changes} bonds are formed/broken when going from start to end.\n") # Determine bond orders at start and end all_bos = [bond_orders_for_geom(geom, bond_indices) for geom in geoms] # Reference bond lengths r0s = get_r0s(start, bond_indices) # Print summary for title, geom, bos in zip(("Start", "End"), geoms, all_bos): print(highlight_text(title, level=1)) for bo_, (a, b) in zip(bos, bond_indices): aa = atoms[a] ba = atoms[b] print(f"Bond ({aa}{a: >4},{ba}{b: >4}): BO={bo_:.2f}") print( f"mean(BOs)={bos.mean():.2f}, max(BOs)={bos.max():.2f}, min(BOs)={bos.min():.2f}" ) print() # Guess goal bond lenghts from mean bond order start_bos, end_bos = all_bos mean_bos = (start_bos + end_bos) / 2 # Determine goal bond lenghts that yield the mean bond orders goal_lengths = length_for_bond_order(mean_bos, r0s) # Constrain formed/broken bonds constrain_prims = normalize_prim_inputs([["BOND", *bond] for bond in bond_indices]) # Start relaxed scans towards TS from 'start' and 'end' print(highlight_text("Relaxed scan towards TS, from 'start'")) start_guess, *_ = relaxed_scan( start, calc_getter, constrain_prims, target_values=goal_lengths, title="start" ) print() start_guess.dump_xyz("start_guess") print(highlight_text("Relaxed scan towards TS, from 'end'")) end_guess, *_ = relaxed_scan( end, calc_getter, constrain_prims, target_values=goal_lengths, title="end" ) print() end_guess.dump_xyz("end_guess") # Compare energies and take the one with lower energy start_en = end_en = print(f"Energy of 'start'-guess: {start_en:.6f} au") print(f"Energy of 'end'-guess: {end_en:.6f} au") start_lower = start_en <= end_en lower = "'start'-guess" if start_lower else "'end'-guess" print(f"{lower} has lower energy.") ts_guess = start_guess if start_lower else end_guess if recreate: print("Recreated Geometry object.") ts_guess = Geometry(ts_guess.atoms, ts_guess.cart_coords, coord_type="redund") ts_guess.set_calculator(calc_getter()) print() print(f"Coordinates of TS guess\n{ts_guess.as_xyz()})\n") ts_guess_fn = "" ts_guess.dump_xyz(ts_guess_fn) print(f"Dumped TS-guess to '{ts_guess_fn}'.") return ts_guess