Source code for pysisyphus.intcoords.update

import numpy as np

from pysisyphus.config import BEND_MIN_DEG, LB_MIN_DEG
from pysisyphus.helpers_pure import log
from pysisyphus.intcoords.eval import eval_primitives
from pysisyphus.intcoords.exceptions import NeedNewInternalsException
from pysisyphus.intcoords.valid import bend_valid, dihedral_valid
from pysisyphus.intcoords.PrimTypes import Bends, Dihedrals, Rotations

[docs] def correct_dihedrals(new_dihedrals, old_dihedrals): """Dihedrals are periodic. Going from -179° to 179° is not a step of 358°, but a step of 2°. By considering the actual distance of the dihedrals from π the correct step can be calculated. dihedral step length = abs(abs(new_dihedral) - π) + abs(abs(old_dihedral) - π) or put differently dihedral step length = abs(abs(new_dihedral - old_dihedral) - 2*π) The sign is left to be determined. Going from -179° to 179° (roughly π - -π = 2π) is a counter clockwise rotation and the dihedral has to decrease below -π. Going from 179° to -179° (roughly -π - π = -2π) is a clockwise rotation and the dihedral increases abvove π. So the correct sign corresponds to the negative sign of the original difference. original difference 2π -> dihedral must decrease -> sign = -1 original difference -2π -> dihedral must increase -> sign = +1 Overall, the old dihedral is modified by the actual step length with the correct sign.""" new_dihedrals = np.atleast_1d(new_dihedrals) old_dihedrals = np.atleast_1d(old_dihedrals) dihedrals_step = new_dihedrals - old_dihedrals shifted_by_2pi = np.abs(np.abs(dihedrals_step) - 2 * np.pi) < np.pi / 2 corrected_dihedrals = new_dihedrals.copy() corrected_dihedrals[shifted_by_2pi] -= ( 2 * np.pi * np.sign(dihedrals_step[shifted_by_2pi]) ) return corrected_dihedrals
[docs] def update_internals( new_coords3d, old_internals, primitives, dihedral_inds, rotation_inds, bend_inds, check_dihedrals=False, check_bends=False, bend_min_deg=BEND_MIN_DEG, bend_max_deg=LB_MIN_DEG, rotation_thresh=0.9, logger=None, ): prim_internals = eval_primitives(new_coords3d, primitives) new_internals = np.array([prim_int.val for prim_int in prim_internals]) internal_diffs = new_internals - old_internals new_rotations = new_internals[rotation_inds] # Check for approaching singularity as discussed in the geomeTRIC paper. The # original code seems to check this only for linear molecules and instead # does some +2π/-2π magic, similar to how dihedrals differences are handled. if (np.abs(new_rotations / np.pi) >= rotation_thresh).any(): raise NeedNewInternalsException(new_coords3d) dihedrals = [prim_internals[i] for i in dihedral_inds] dihedral_diffs = internal_diffs[dihedral_inds] # Find differences that are shifted by 2*pi shifted_by_2pi = np.abs(np.abs(dihedral_diffs) - 2 * np.pi) < np.pi / 2 new_dihedrals = np.array([dihed.val for dihed in dihedrals]) if any(shifted_by_2pi): new_dihedrals[shifted_by_2pi] -= ( 2 * np.pi * np.sign(dihedral_diffs[shifted_by_2pi]) ) # Update values for dihed, new_val in zip(dihedrals, new_dihedrals): dihed.val = new_val invalid_inds = list() # See if dihedrals became invalid (collinear atoms) if check_dihedrals: are_valid = [dihedral_valid(new_coords3d, prim.inds) for prim in dihedrals] try: first_dihedral = dihedral_inds[0] except IndexError: first_dihedral = 0 invalid_inds = [ i + first_dihedral for i, is_valid in enumerate(are_valid) if not is_valid ] if check_bends and len(bend_inds) > 0: bends = [prim_internals[i] for i in bend_inds] are_valid = [ bend_valid(new_coords3d, prim.inds, bend_min_deg, bend_max_deg) for prim in bends ] first_bend = bend_inds[0] invalid_inds = [ i + first_bend for i, is_valid in enumerate(are_valid) if not is_valid ] if len(invalid_inds) > 0: invalid_prims = [primitives[i] for i in invalid_inds] invalid_msg = ", ".join([str(tp) for tp in invalid_prims]) log(logger, "Internal coordinate(s) became invalid! Need new internal coordinates!") log(logger, f"Invalid primitives: {invalid_msg}") raise NeedNewInternalsException( new_coords3d, invalid_inds=invalid_inds, invalid_prims=invalid_prims ) return prim_internals
[docs] def inds_from_prim_types(typed_prims, prim_types): inds = [ i for i, (prim_type, *inds) in enumerate(typed_prims) if prim_type in prim_types ] return inds
[docs] def transform_int_step( int_step, old_cart_coords, cur_internals, Bt_inv_prim, primitives, typed_prims=None, dihedral_inds=None, rotation_inds=None, bend_inds=None, check_dihedrals=False, check_bends=False, bend_min_deg=BEND_MIN_DEG, bend_max_deg=LB_MIN_DEG, freeze_atoms=None, constrained_inds=None, update_constraints=False, cart_rms_thresh=1e-6, Bt_inv_prim_getter=None, max_cycles=25, logger=None, ): """Transformation is done in primitive internals, so int_step must be given in primitive internals and not in DLC!""" # If an iterable of typed prims is given we can derive bend/dihedral/rotation # indices from them. if typed_prims is not None: if dihedral_inds is None: dihedral_inds = inds_from_prim_types(typed_prims, Dihedrals) if bend_inds is None: bend_inds = inds_from_prim_types(typed_prims, Bends) if rotation_inds is None: rotation_inds = inds_from_prim_types(typed_prims, Rotations) if freeze_atoms is None: freeze_atoms = list() if constrained_inds is None: constrained_inds = list() freeze_atoms = np.array(freeze_atoms, dtype=int) new_cart_coords = old_cart_coords.copy() remaining_int_step = int_step target_internals = cur_internals + int_step # When we want to update the constraints we use the target primitive internals, if update_constraints: constrained_vals = target_internals[constrained_inds] # otherwise we try to stay with the original constrained values. else: constrained_vals = cur_internals[constrained_inds] def backtransform(remaining_int_step, Bt_inv_prim): """Separate function so it can be a called after the main loop finished to re-enforce the constraints.""" nonlocal new_cart_coords cart_step = # Remove step from frozen atoms. cart_step.reshape(-1, 3)[freeze_atoms] = 0.0 cart_rms = np.sqrt(np.mean(cart_step**2)) # Update cartesian coordinates new_cart_coords += cart_step # Determine new internal coordinates new_prim_ints = update_internals( new_cart_coords.reshape(-1, 3), old_internals, primitives, dihedral_inds, rotation_inds, bend_inds, check_dihedrals=check_dihedrals, check_bends=check_bends, bend_min_deg=bend_min_deg, bend_max_deg=bend_max_deg, logger=logger, ) new_internals = [prim.val for prim in new_prim_ints] return new_prim_ints, new_internals, cart_rms last_rms = 9999 old_internals = cur_internals backtransform_failed = True for i in range(max_cycles): if Bt_inv_prim_getter is not None: Bt_inv_prim = Bt_inv_prim_getter(new_cart_coords) log(logger, f"Recalculated (B^T)^+ in microcycle {i}") new_prim_ints, new_internals, cart_rms = backtransform( remaining_int_step, Bt_inv_prim ) remaining_int_step = target_internals - new_internals internal_rms = np.sqrt(np.mean(remaining_int_step**2)) log( logger, f"Cycle {i}: rms(Δcart)={cart_rms:1.4e}, rms(Δint.) = {internal_rms:1.5e}", ) # This assumes the first cart_rms won't be > 9999 ;) if cart_rms < last_rms: # Store results of the conversion cycle for laster use, if # the internal-cartesian-transformation goes bad. best_cycle = (new_cart_coords.copy(), new_internals.copy()) best_cycle_ind = i elif i != 0: # If the conversion somehow fails we fallback to the best previous step. log(logger, f"Backconversion failed! Falling back to step {best_cycle_ind}") new_cart_coords, new_internals = best_cycle break else: raise Exception( "Internal-cartesian back-transformation already " "failed in the first step. Aborting!" ) old_internals = new_internals last_rms = cart_rms if cart_rms < cart_rms_thresh: log( logger, f"Internal->Cartesian transformation converged in {i+1} cycle(s)!", ) backtransform_failed = False break if len(constrained_inds) > 0: for j in range(max_cycles): cur_constrained_vals = np.array(new_internals)[constrained_inds] diff = constrained_vals - cur_constrained_vals if any(np.abs(diff) <= 1e-5): break remaining_int_step = np.zeros_like(remaining_int_step) remaining_int_step[constrained_inds] = diff new_prim_ints, new_internals, _ = backtransform( remaining_int_step, Bt_inv_prim ) if j > 0: log(logger, f"Re-enforced constraints in {j} additional cycle(s).") log(logger, "") # Return the difference between the new cartesian coordinates that yield # the desired internal coordinates and the old cartesian coordinates. cart_step = new_cart_coords - old_cart_coords return new_prim_ints, cart_step, backtransform_failed