tscode - pypi Package Compare versions

Comparing version 0.4.9 to 0.4.11

CHANGELOG.md

@@ -8,2 +8,12 @@ ## TSCoDe Changelog
 
+## 0.4.11 (March 10, 2024)
+- Various small bugfixes/graphic restyling
+- Stuctures are optimized at embedder.options.theory_level before running the mtd_search> operator, and a FatalError is raised if they scramble during this process. This increases the robustness of the workflow by avoiding changes in the molecular graph if the input structure is bad.
+- Added a warning printout system to the Embedder class, via logging and appending strings to the embedder.warnings list.
+- Added a saturation and compenetration check for input molecules, that warns about potential bad geometries. 
+
+## 0.4.10 (February 24, 2024)
+- Reduced/summarized printouts for loose to tight optimization steps. 
+<!-- - Removed call to compenetration_refining() on input ensemble: it is usually already pruned for compenetration for embed runs, and other ensembles (for example from "refine> mtd>" runs coming from CREST) benefit from relaxing eventual clashes that are present. -->
+
 ## 0.4.9 (February 22, 2024)
@@ -10,0 +20,0 @@ - CREST constraints (mtd_search>) are now passed as distance constraints instead of atom fixing (forgot to implement it this way before?).

docs/examples.rst

@@ -50,3 +50,3 @@ .. _exs:
 
-    DIST(a=2.135)
+    DIST(A=2.135)
 
@@ -53,0 +53,0 @@ maleimide.xyz 0A 5x

PKG-INFO

 Metadata-Version: 2.1
 Name: tscode
-Version: 0.4.9
+Version: 0.4.11
 Summary: Computational chemistry general purpose transition state builder and ensemble optimizer
@@ -43,2 +43,12 @@ Home-page: https://www.github.com/ntampellini/tscode
 
+## 0.4.11 (March 10, 2024)
+- Various small bugfixes/graphic restyling
+- Stuctures are optimized at embedder.options.theory_level before running the mtd_search> operator, and a FatalError is raised if they scramble during this process. This increases the robustness of the workflow by avoiding changes in the molecular graph if the input structure is bad.
+- Added a warning printout system to the Embedder class, via logging and appending strings to the embedder.warnings list.
+- Added a saturation and compenetration check for input molecules, that warns about potential bad geometries. 
+
+## 0.4.10 (February 24, 2024)
+- Reduced/summarized printouts for loose to tight optimization steps. 
+<!-- - Removed call to compenetration_refining() on input ensemble: it is usually already pruned for compenetration for embed runs, and other ensembles (for example from "refine> mtd>" runs coming from CREST) benefit from relaxing eventual clashes that are present. -->
+
 ## 0.4.9 (February 22, 2024)
@@ -45,0 +55,0 @@ - CREST constraints (mtd_search>) are now passed as distance constraints instead of atom fixing (forgot to implement it this way before?).

tscode.egg-info/PKG-INFO

 Metadata-Version: 2.1
 Name: tscode
-Version: 0.4.9
+Version: 0.4.11
 Summary: Computational chemistry general purpose transition state builder and ensemble optimizer
@@ -43,2 +43,12 @@ Home-page: https://www.github.com/ntampellini/tscode
 
+## 0.4.11 (March 10, 2024)
+- Various small bugfixes/graphic restyling
+- Stuctures are optimized at embedder.options.theory_level before running the mtd_search> operator, and a FatalError is raised if they scramble during this process. This increases the robustness of the workflow by avoiding changes in the molecular graph if the input structure is bad.
+- Added a warning printout system to the Embedder class, via logging and appending strings to the embedder.warnings list.
+- Added a saturation and compenetration check for input molecules, that warns about potential bad geometries. 
+
+## 0.4.10 (February 24, 2024)
+- Reduced/summarized printouts for loose to tight optimization steps. 
+<!-- - Removed call to compenetration_refining() on input ensemble: it is usually already pruned for compenetration for embed runs, and other ensembles (for example from "refine> mtd>" runs coming from CREST) benefit from relaxing eventual clashes that are present. -->
+
 ## 0.4.9 (February 22, 2024)
@@ -45,0 +55,0 @@ - CREST constraints (mtd_search>) are now passed as distance constraints instead of atom fixing (forgot to implement it this way before?).

tscode/__main__.py

@@ -26,3 +26,3 @@ # coding=utf-8
 
-__version__ = '0.4.9'
+__version__ = '0.4.11'
 
@@ -29,0 +29,0 @@ if __name__ == '__main__':

tscode/automep.py

 import time
 
+from ase import Atoms
 import numpy as np
@@ -9,2 +10,4 @@ from networkx import cycle_basis
 from tscode.graph_manipulations import neighbors
+from tscode.hypermolecule_class import align_structures
+from tscode.mep_relaxer import interpolate_structures
 from tscode.optimization_methods import optimize
@@ -14,5 +17,6 @@ from tscode.utils import graphize, write_xyz
 
-def automep(embedder, n_images=7):
+def automep(embedder, n_images=9):
 
     assert embedder.options.calculator == "XTB"
+    assert len(embedder.objects) == 2, "Provide two molecules as start/endpoints."
 
@@ -22,8 +26,8 @@ mol = embedder.objects[0]
 
-    # Get cycle indices bigger than 6
+    # Get cycle indices between 7 and 9
     graph = graphize(coords, mol.atomnos)
-    cycles = [l for l in cycle_basis(graph) if len(l) == 7]
-    assert len(cycles) == 1, "Automep only works for 7-membered ring flips at the moment"
+    cycles = [l for l in cycle_basis(graph) if len(l) in (7, 8, 9)]
+    assert len(cycles) == 1, "Automep only works for 7/8/9-membered ring flips at the moment"
 
-    embedder.log('--> AutoMEP - Building MEP for 7-membered ring inversion')
+    embedder.log(f'--> AutoMEP - Building MEP for {len(cycles[0])}-membered ring inversion')
     embedder.log(f'    Preoptimizing starting point at {embedder.options.calculator}/{embedder.options.theory_level}({embedder.options.solvent}) level')
@@ -45,24 +49,55 @@
 
-    start_angles = np.array([dihedral(coords[d]) for d in dihedrals+exocyclic])
+    # start_angles = np.array([dihedral(coords[d]) for d in dihedrals+exocyclic])
     target_angles = np.array([0 for _ in dihedrals] + [180 for _ in exocyclic])
-    multipliers = np.linspace(1, -1, n_images)
+    print("Optimizing planar TS guess...", end="\r")
+    ts_guess, _, _ = xtb_opt(
+                            coords,
+                            mol.atomnos,
+                            constrained_dihedrals=dihedrals+exocyclic,
+                            constrained_dih_angles=target_angles,
+                            method=embedder.options.theory_level,
+                            solvent=embedder.options.solvent,
+                            procs=embedder.procs
+                        )
+    # multipliers = np.linspace(1, -1, n_images)
 
-    mep_angles = [(start_angles * m + target_angles * (1-m)) % 360 for m in multipliers]
+    # mep_angles = [(start_angles * m + target_angles * (1-m)) % 360 for m in multipliers]
+    # mep = []
+    # for i, m_a in enumerate(mep_angles):
+    #     t_start = time.perf_counter()
+    #     coords, _, _ = xtb_opt(coords,
+    #                             mol.atomnos,
+    #                             constrained_dihedrals=dihedrals+exocyclic,
+    #                             constrained_dih_angles=m_a,
+    #                             method=embedder.options.theory_level,
+    #                             solvent=embedder.options.solvent,
+    #                             procs=embedder.procs)
+    #     embedder.log(f'    - optimized image {i+1}/{len(mep_angles)} ({round(time.perf_counter()-t_start, 3)} s)')
+    #     mep.append(coords)
 
-    mep = []
-    for i, m_a in enumerate(mep_angles):
-        t_start = time.perf_counter()
-        coords, _, _ = xtb_opt(coords,
+    mep = interpolate_structures(align_structures(np.array([coords,
+                                                            ts_guess,
+                                                            embedder.objects[1].atomcoords[0]])),
+                                 mol.atomnos,
+                                 n=n_images,
+                                 method='linear')
+    
+    for g, geom in enumerate(mep):
+        if g not in (0, n_images-1):
+            print(f"Relaxing image {g+1}/{n_images}...", end="\r")
+            positions = geom.get_positions()
+            opt_geom, _, _ = xtb_opt(
+                                positions,
                                 mol.atomnos,
                                 constrained_dihedrals=dihedrals+exocyclic,
-                                constrained_dih_angles=m_a,
+                                constrained_dih_angles=[dihedral(positions[quadruplet]) for quadruplet in dihedrals+exocyclic],
                                 method=embedder.options.theory_level,
                                 solvent=embedder.options.solvent,
-                                procs=embedder.procs)
-        embedder.log(f'    - optimized image {i+1}/{len(mep_angles)} ({round(time.perf_counter()-t_start, 3)} s)')
-        mep.append(coords)
+                                procs=embedder.procs
+                            )
+            mep[g] = Atoms(mol.atomnos, positions=opt_geom)
 
     with open(f"{mol.rootname}_automep.xyz", "w") as f:
         for c in mep:
-            write_xyz(c, mol.atomnos, f)
+            write_xyz(c.get_positions(), mol.atomnos, f)
 
@@ -69,0 +104,0 @@ embedder.log(f"\n--> Saved autogenerated MEP as {mol.rootname}_automep.xyz\n")

tscode/embedder_options.py

@@ -304,3 +304,3 @@ # coding=utf-8
         if self.keywords_simple:
-            embedder.log('--> Parsed keywords, in order of execution:\n    ' + ' '.join(self.sorted_keywords()) + '\n')
+            embedder.log('\n--> Parsed keywords, in order of execution:\n    ' + ' '.join(self.sorted_keywords()) + '\n')
 
@@ -307,0 +307,0 @@ def refine(self, options, *args):

tscode/errors.py

@@ -55,2 +55,7 @@ # coding=utf-8
     Thrown when trying to access orbital data when they are not present
+    '''
+
+class FatalError(Exception):
+    '''
+    Thrown when a molecule optimization crashed or scrambled fatally
     '''

tscode/operators.py

@@ -31,3 +31,3 @@ # coding=utf-8
 from tscode.calculators._xtb import crest_mtd_search
-from tscode.errors import InputError
+from tscode.errors import InputError, FatalError
 from tscode.graph_manipulations import graphize
@@ -45,3 +45,3 @@ from tscode.hypermolecule_class import align_structures
                                    prune_conformers_rmsd_rot_corr)
-from tscode.utils import (get_scan_peak_index, read_xyz,
+from tscode.utils import (get_scan_peak_index, molecule_check, read_xyz,
                           suppress_stdout_stderr, time_to_string, write_xyz)
@@ -452,3 +452,46 @@
     logfunction = embedder.log
+    constrained_indices = _get_internal_constraints(filename, embedder)
+    constrained_distances = [embedder.get_pairing_dists_from_constrained_indices(cp) for cp in constrained_indices]
 
+    logfunction(f'--> {filename}: Geometry optimization pre-mtd_search ({embedder.options.theory_level} via {embedder.options.calculator})')
+    logfunction(f'    {len(constrained_indices)} constraints applied{": "+constrained_indices if len(constrained_indices) > 0 else ""}')
+    
+    for c, coords in enumerate(mol.atomcoords.copy()):
+        logfunction(f"    Optimizing conformer {c+1}/{len(mol.atomcoords)}")
+
+        opt_coords, _, success = optimize(
+                                    coords,
+                                    mol.atomnos,
+                                    calculator=embedder.options.calculator,
+                                    method=embedder.options.theory_level,
+                                    solvent=embedder.options.solvent,
+                                    charge=embedder.options.charge,
+                                    procs=embedder.procs,
+                                    constrained_indices=constrained_indices,
+                                    constrained_distances=constrained_distances,
+                                    title=f'{filename.split(".")[0]}_conf{c+1}',
+                                ) if embedder.options.optimization else coords
+        
+        exit_status = "" if success else "CRASHED"
+        
+        if success:
+            success = molecule_check(coords, opt_coords, mol.atomnos)
+            exit_status = "" if success else "SCRAMBLED"
+
+        if not success:
+            dumpname = filename.split(".")[0] + f"_conf{c+1}_{exit_status}.xyz"
+            with open(dumpname, "w") as f:
+                write_xyz(opt_coords, mol.atomnos, f, title=f"{filename}, conformer {c+1}/{len(mol.atomcoords)}, {exit_status}")
+
+            logfunction(f"{filename}, conformer {c+1}/{len(mol.atomcoords)} optimization {exit_status}. Inspect geometry at {dumpname}. Aborting run.")
+
+            raise FatalError(filename)
+        
+        # update embedder structures after optimization
+        mol.atomcoords[c] = opt_coords
+
+    # update mol and embedder graph after optimization 
+    mol.graph = graphize(mol.atomcoords[0], mol.atomnos)
+    embedder.graphs = [m.graph for m in embedder.objects]
+
     max_workers = embedder.avail_cpus//2 or 1
@@ -472,4 +515,2 @@ logfunction(f'--> Performing {embedder.options.calculator} GFN2//GFN-FF' + (
         t_start_conf = time.perf_counter()
-        constrained_indices = _get_internal_constraints(filename, embedder)
-        constrained_distances = [embedder.get_pairing_dists_from_constrained_indices(cp) for cp in constrained_indices]
         try:
@@ -476,0 +517,0 @@ conf_batch = crest_mtd_search(

tscode/solvents.py

@@ -89,2 +89,3 @@ # coding=utf-8
     'dcm':'ch2cl2',
+    'dichloromethane':'ch2cl2',
     'carbondisuphide':'cs2',
@@ -91,0 +92,0 @@ 'carbondisulfide':'cs2',

tscode/torsion_module.py

@@ -1156,4 +1156,7 @@ # coding=utf-8
     for hb in hydrogen_bonds:
-        graph.remove_edge(*hb)
+        try:
+            graph.remove_edge(*hb)
+        except nx.NetworkXError():
+            pass
 
     return structures[final_mask], final_mask

tscode/utils.py

@@ -355,3 +355,3 @@ # coding=utf-8
 
-def scramble_check(TS_structure, TS_atomnos, excluded_atoms, mols_graphs, max_newbonds=0) -> bool:
+def scramble_check(TS_structure, TS_atomnos, excluded_atoms, mols_graphs, max_newbonds=0, logfunction=False, title=None) -> bool:
     '''
@@ -384,2 +384,4 @@ Check if a multimolecular arrangement has scrambled during some optimization
     if len(delta_bonds) > max_newbonds:
+        if logfunction is not None:
+            logfunction(f"{title}, scramble_check - found {len(delta_bonds)} extra bonds: {delta_bonds}")
         return False
@@ -467,2 +469,43 @@
     
-    return func_return, payload, elapsed
+    return func_return, payload, elapsed
+
+def _saturation_check(atomnos, charge=0, embedder=None):
+
+    transition_metals = [
+                    "Sc", "Ti", "V", "Cr", "Mn", "Fe",
+                    "Co", "Ni", "Cu", "Zn", "Y", "Zr",
+                    "Nb", "Mo", "Tc", "Ru", "Rh", "Pd",
+                    "Ag", "Cd", "La", "Ce", "Pr", "Nd",
+                    "Pm", "Sm", "Eu", "Gd", "Tb", "Dy",
+                    "Ho", "Er", "Tm", "Yb", "Lu", "Hf",
+                    "Ta", "W", "Re", "Os", "Ir", "Pt",
+                    "Au", "Hg", "Th", "Pa", "U", "Np",
+                    "Pu", "Am", 
+    ]
+
+    # if we have any transition metal, it's hard to tell
+    # if the structure looks ok: in this case we assume it is.
+    organometallic = any([pt[a].symbol in transition_metals for a in atomnos])
+
+    odd_valent = [  #1 valent
+                    "H", "Li", "Na", "K", "Rb", "Cs",
+                    "F", "Cl", "Br", "I", "At",
+
+                    # 3/5 valent
+                    "N", "P", "As", "Sb", "Bi",
+                    "B", "Al", "Ga", "In", "Tl",
+                 ]
+
+    n_odd_valent = sum([1 for a in atomnos if pt[a].symbol in odd_valent])
+    looks_ok = ((n_odd_valent + charge) / 2) % 1 < 0.001 if not organometallic else True
+
+    if embedder is not None:
+        if looks_ok:
+            embedder.log(f"--> Saturation check passed: input structure has an even saturation index")
+
+        else:
+            s = f"--> WARNING! Saturation check failed: input structure has an odd saturation index (charge={charge}). Radical or bad input geometry?"
+            embedder.log(s)
+            embedder.warnings.append(s)
+
+    return looks_ok

Alert delta unavailable

Currently unable to show alert delta for PyPI packages.

Improved metrics

Total package byte prevSize

1362936

1.22%

Lines of code

18933

0.8%