20 adding vibration mode analysis

TCutility/__init__.py

@@ -1,4 +1,4 @@
 ensure_list = lambda x: [x] if not isinstance(x, (list, tuple, set)) else x  # noqa: E731
 squeeze_list = lambda x: x[0] if len(x) == 1 else x  # noqa: E731
 
-from TCutility import results, constants  # noqa: F401, E402
+from TCutility import analysis, results, constants  # noqa: F401, E402

TCutility/analysis/ts_vibration.py

@@ -0,0 +1,139 @@
+import numpy as np
+import warnings
+from TCutility import results
+from scm import plams
+
+def avg_relative_bond_length_delta(base: plams.Molecule, pos: plams.Molecule, neg: plams.Molecule, atom1: int, atom2: int) -> float:
+    '''Function to calculate relative atom distance change in vibrational mode
+
+    Args:
+        base: plams.Molecule object containing the coordinates before applying vibration
+        pos: plams.Molecule object with the vibrational mode added
+        neg: plams.Molecule object with the vibrational mode subtracted
+        atom1: label for the first atom
+        atom2: label for the second atom
+
+    Returns:
+        Average relative difference from the baseline distance between selected atoms in this vibrational mode (as percentage). 
+    ''' 
+    basedist = base[atom1].distance_to(base[atom2])
+    x = abs((pos[atom1].distance_to(pos[atom2])/basedist)-1) # distances between positive and negative vibration are normalized with base distance
+    y = abs((neg[atom1].distance_to(neg[atom2])/basedist)-1) # variances are obtained and averaged
+    return (x+y)/2
+
+def determine_ts_reactioncoordinate(data: results.Result, mode_index: int = 0, bond_tolerance: float = 1.28, min_delta_dist: float = 0.0) -> np.ndarray:
+    '''Function to retrieve reaction coordinate from a given transitionstate, using the first imaginary frequency.
+
+    Args:
+        data: TCutility.results.Result object containing calculation data
+        mode_index: vibrational mode index to analyze
+        bond_tolerance: parameter for plams.Molecule.guess_bonds() function
+        min_delta_dist: minimum relative bond length change before qualifying as active atom. If 0, all bond changes are counted
+
+    Returns:
+        Array containing all the obtained reaction coordinates. Reaction coordinate format is [active_atom1, active_atom2, sign], using Distance reactioncoordinate
+        Symmetry elements are ignored, by convention the atom labels are increasing (atom1 < atom2)
+    '''
+
+    assert 'modes' in data.properties.vibrations, 'Vibrational data is required, but was not present in .rkf file'
+
+    outputmol = data.molecule.output
+
+    base = np.array(outputmol).reshape(-1,3)
+    tsimode = np.array(data.properties.vibrations.modes[mode_index]).reshape(-1,3)
+
+    posvib = outputmol.copy()
+    posvib.from_array(base + tsimode)
+    posvib.guess_bonds(dmax=bond_tolerance)
+    negvib = outputmol.copy()
+    negvib.from_array(base - tsimode)
+    negvib.guess_bonds(dmax=bond_tolerance)
+
+    pairs = np.where(posvib.bond_matrix() - negvib.bond_matrix()) 
+    rc = np.array(
+                [   [a+1, b+1, np.sign(posvib.bond_matrix()[a][b] - negvib.bond_matrix()[a][b])]
+                    for a, b in np.c_[pairs[0], pairs[1]]  # bond matrices and pairs matrix are 0-indexed unlike the plams.Molecule labels, hence a+1 and b+1 are needed
+                    if a < b and avg_relative_bond_length_delta(outputmol, posvib, negvib, a+1, b+1) > min_delta_dist
+                ],  dtype=int)
+    return rc
+
+
+def validate_transitionstate(calc_dir: str, rcatoms: list = None, analyze_modes: int = 1, **kwargs) -> bool:
+    ''' Function to determine whether a transition state calculation yielded the expected transition state. Checks the reaction coordinates provided by the user (or in the .rkf file User Input section) and compares
+        this against the reaction coordinates found in the imaginary modes of the transitionstate. If the transitionstate has multiple imaginary frequencies, it is possible to check multiple modes for the expected
+        reaction coordinate.
+
+    Args:
+        calc_dir: path pointing to the desired calculation.
+        rcatoms: list or array containing expected reaction coordinates, to check against the transition state. If not defined, it is obtained from the ams.rkf user input. 
+                  Only uses 'Distance' reaction coordinate. Format should be [atomlabel1, atomlabel2, (optional) sign]
+        analyze_modes: Number of imaginary modes to analyze, default only the first mode. Modes are ordered lowest frequency first. If 0 or negative value is provided, analyze all modes with imaginary frequency.
+        **kwargs: keyword arguments for use in :func:`determine_ts_reactioncoordinate`.
+
+    Returns:
+        Boolean value, True if the found transition state reaction coordinates contain the expected reaction coordinates, False otherwise. 
+        If multiple modes are analyzed, returns True if at least one mode contains the expected reaction coordinates.
+    '''
+
+    data = results.read(calc_dir)
+    assert 'modes' in data.properties.vibrations, 'Vibrational data is required, but was not present in .rkf file'
+
+    nimag = sum(1 for x in data.properties.vibrations.frequencies if x < 0)
+    if nimag == 0: 
+        return False # no imaginary modes found in transitionstate
+
+    if analyze_modes > 0:                   # if positive value is provided, check that many imaginary vibrational modes. Otherwise all imaginary ones are analyzed
+        nimag = min(nimag, analyze_modes)   # limit to the actual number of vibrational modes
+
+    if isinstance(rcatoms, type(None)):
+        assert 'reactioncoordinate' in data.input.transitionstatesearch, 'Reaction coordinate is a required input, but was neither provided nor present in the .rkf file'
+        rcatoms = np.array([[int(float(y)) for y in x.split()[1:]] for x in data.input.transitionstatesearch.reactioncoordinate if x.split()[0] == 'Distance'])
+        assert len(rcatoms) > 0, 'Reaction coordinate data was present in .rkf file, but no reaction coordinate using Distance was provided'
+
+    # cast the reaction coordinate as 2d array. If reaction coordinate sign is provided, it must be provided for all reaction coordinates. 
+    # the np.atleast_2d() function will cause a deprecation warning if the dimensions mismatch, this is raised as a value error instead
+    with warnings.catch_warnings(record=True) as w:
+        rcatoms = np.atleast_2d(rcatoms)
+        if len(w) > 0:
+            raise ValueError(w[-1].message)
+
+    assert np.all([len(x)==2 or len(x)==3 for x in rcatoms]), 'Invalid format of reaction coordinate. Reaction coordinate format must be [label1, label2, (optional) sign]'
+
+    ret = []
+
+    for idx in range(nimag):
+        result = determine_ts_reactioncoordinate(data, mode_index=idx, **kwargs)
+        if len(result) == 0:
+            ret.append(False)
+            continue
+
+        # sort for consistency
+        for index, [a,b,*c] in enumerate(rcatoms):
+            if c:
+                c = np.sign(c)
+            if a > b:
+                a,b = b,a
+            rcatoms[index] = [a,b,*c]
+        rcatoms = rcatoms[rcatoms[:,1].argsort()]
+        rcatoms = rcatoms[rcatoms[:,0].argsort()]
+        result = result[result[:,1].argsort()]
+        result = result[result[:,0].argsort()]
+
+        # (optional) check for internal consistency if the reaction coordinate sign is provided
+        # only the first reaction coordinate sign is arbitrary, check remaining coordinates for consistency
+        if len(rcatoms[0]) > 2:
+            foundmatch = False
+            for idx, rc in enumerate(result):
+                if rc[0] == rcatoms[0][0] and rc[1] == rcatoms[0][1]:
+                    result[:, 2] = result[:, 2] if np.sign(result[idx][2]) == rcatoms[0][2] else -result[:, 2]
+                    foundmatch = True
+                    break
+            if not foundmatch: # at least one element of rcatoms is not present in result
+                ret.append(False)
+                continue
+        else:
+            result = result[:, :-1]
+
+        ret.append(set(str(x) for x in rcatoms).issubset(set(str(y) for y in result)))
+
+    return any(ret)

TCutility/constants.py

@@ -1 +1,3 @@
 HA2KCALMOL: float = 627.509474 
+ANG2BOHR: float = 1.8897259886
+BOHR2ANG: float = 0.5291772490

TCutility/results/adf.py

@@ -45,12 +45,11 @@ def get_calc_settings(info: Result) -> Result:
     # determine if SFOs are unrestricted or not
     ret.unrestricted_sfos = ('SFOs', 'energy_B') in reader_adf
 
-    if ('Geometry', 'grouplabel') in reader_adf:
-        ret.symmetry.group = reader_adf.read('Geometry', 'grouplabel').strip()
+    # get the symmetry group
+    ret.symmetry.group = reader_adf.read('Geometry', 'grouplabel').strip()
 
     # get the symmetry labels
-    if ('Symmetry', 'symlab') in reader_adf:
-        ret.symmetry.labels = reader_adf.read('Symmetry', 'symlab').strip().split()
+    ret.symmetry.labels = reader_adf.read('Symmetry', 'symlab').strip().split()
 
     # determine if MOs are unrestricted or not
     ret.unrestricted_mos = (ret.symmetry.labels[0], 'eps_B') in reader_adf
@@ -127,11 +126,12 @@ def read_vibrations(reader: cache.TrackKFReader) -> Result:
         ret.vibrations = read_vibrations(reader_adf)
 
     # read the Voronoi Deformation Charges Deformation (VDD) before and after SCF convergence (being "inital" and "SCF")
-    vdd_scf: List[float] = ensure_list(reader_adf.read('Properties', 'AtomCharge_SCF Voronoi'))  # type: ignore since plams does not include typing for KFReader. List[float] is returned
-    vdd_ini: List[float] = ensure_list(reader_adf.read('Properties', 'AtomCharge_initial Voronoi'))  # type: ignore since plams does not include typing for KFReader. List[float] is returned
+    if 'Properties' in reader_adf:
+        vdd_scf: List[float] = ensure_list(reader_adf.read('Properties', 'AtomCharge_SCF Voronoi'))  # type: ignore since plams does not include typing for KFReader. List[float] is returned
+        vdd_ini: List[float] = ensure_list(reader_adf.read('Properties', 'AtomCharge_initial Voronoi'))  # type: ignore since plams does not include typing for KFReader. List[float] is returned
 
-    # VDD charges are scf - initial charges. Note, these are in units of electrons while most often these are denoted in mili-electrons
-    ret.vdd.charges = [float((scf - ini)) for scf, ini in zip(vdd_scf, vdd_ini)]
+        # VDD charges are scf - initial charges. Note, these are in units of electrons while most often these are denoted in mili-electrons
+        ret.vdd.charges = [float((scf - ini)) for scf, ini in zip(vdd_scf, vdd_ini)]
 
     # Possible enhancement: get the VDD charges per irrep, denoted by the "Voronoi chrg per irrep" in the "Properties" section in the adf.rkf. 
     # The ordering is not very straightfoward so this is a suggestion for the future with keys: ret.vdd.[IRREP]

TCutility/results/ams.py

@@ -1,7 +1,7 @@
 import numpy as np
 from scm import plams
 from TCutility.results import cache, Result
-from TCutility import ensure_list
+from TCutility import constants, ensure_list
 import os
 from datetime import datetime
 import re
@@ -287,7 +287,7 @@ def get_molecules(calc_dir: str) -> Result:
     # read input molecule
     ret.input = plams.Molecule()
     # read in the coordinates, they are given in Bohr, so convert them to Angstrom
-    coords = np.array(reader_ams.read('InputMolecule', 'Coords')).reshape(natoms, 3) * 0.529177
+    coords = np.array(reader_ams.read('InputMolecule', 'Coords')).reshape(natoms, 3) * constants.BOHR2ANG
     # add the atoms to the molecule
     for atnum, coord in zip(atnums, coords):
         ret.input.add_atom(plams.Atom(atnum=atnum, coords=coord))
@@ -301,7 +301,7 @@ def get_molecules(calc_dir: str) -> Result:
 
     # read output molecule
     ret.output = plams.Molecule()
-    coords = np.array(reader_ams.read('Molecule', 'Coords')).reshape(natoms, 3) * 0.529177
+    coords = np.array(reader_ams.read('Molecule', 'Coords')).reshape(natoms, 3) * constants.BOHR2ANG
     for atnum, coord in zip(atnums, coords):
         ret.output.add_atom(plams.Atom(atnum=atnum, coords=coord))
     if ('Molecule', 'fromAtoms') in reader_ams and ('Molecule', 'toAtoms') in reader_ams and ('Molecule', 'bondOrders'):