Breaking bonds in a ring molecule and modifying distances for NEB
*****************************************************************

Building structures Technical DFTB NEB

**Requires:** AMS2026 or later

-[ Related examples ]-

* Nudged Elastic Band (NEB) for gasphase reaction with DFTB

* Automated Transition-State Workflow for a Diels-Alder Addition

* Band structure with DFTB (SCC-DFTB, GFN1-xTB)

-[ Related tutorials ]-

* NEB tutorial

* Building molecules tutorial

-[ Related documentation ]-

* NEB task

* ChemicalSystem geometry manipulation

[image]


Create the systems
==================


Create the ring molecule from SMILES
------------------------------------

First, we create "ring_mol" from the SMILES code, and then we cut the
O-C and N-C bonds and increase the distance between the fragments. We
need to first cut the bonds before increasing the distance, since
"ChemicalSystem.set_distance" requires separate molecule fragments to
work on.

   from scm.plams import view
   from scm.base import ChemicalSystem

   # from_smiles creates bonds between atoms
   ring_mol = ChemicalSystem.from_smiles("C1CON=N1")
   view(
       ring_mol,
       width=350,
       height=350,
       direction="tilt_pca3",
       show_atom_labels=True,
       atom_label_type="Name",
   )

[image]


Identify which bonds to cut
---------------------------

It’s easiest if you know the indices of the bonded atoms. Then you can
just cut those bonds, e.g.
"ring_mol.bonds.remove_bonds_between_atoms(0, 4)". That would remove
the N-C bond in the above picture (there the indices start with 1, but
when calling Python methods they start from 0)

That is fine and the most straightforward way to cut the bonds.
However, it is difficult to automate because the atomic indices may
change for different structures.

Here we show a way to programmatically identify “an N atom that is
bonded to another N and C”, and “an O atom that is bonded to an N and
C”. Then we cut those corresponding bonds.

   from typing import List, Tuple, Iterable, Optional


   def find_first_matching(
       cs: ChemicalSystem, elem: str, neighbors: Iterable[str]
   ) -> Optional[Tuple[int, List[int], List[str]]]:
       """Finds the first atom with element ``elem`` and neighbors ``neighbors``.

       Returns a 3-tuple: (index_of_atom, list_of_neighbor_indices, list_of_neighbor_elements)

       In the return value the list_of_neighbor_indices and list_of_neighbor_elements
       have the same length and correspond to each other.
       """

       sorted_target_neighbors = sorted(neighbors)
       for i, at in enumerate(cs):
           if at.symbol != elem:
               continue
           my_neighbors = list(cs.bonds.get_bonded_atoms(i))
           my_neighbors_elems = [cs.atoms[j].symbol for j in my_neighbors]
           if sorted(my_neighbors_elems) == sorted_target_neighbors:
               return i, my_neighbors, my_neighbors_elems

       return None

   sep_mol = ring_mol.copy()

   # N bonded to N and C
   iN, N_neighs, N_neighs_elems = find_first_matching(sep_mol, "N", ["N", "C"])
   iN_C = N_neighs[N_neighs_elems.index("C")]  # index of neighboring C atom

   # O bonded to N and C
   iO, O_neighs, O_neighs_elems = find_first_matching(sep_mol, "O", ["N", "C"])
   iO_C = O_neighs[O_neighs_elems.index("C")]  # index of neighboring C atom

   print(f"I will cut the bonds between atoms {iN}-{iN_C} and {iO}-{iO_C} (zero-based indexing)")

   I will cut the bonds between atoms 4-0 and 2-1 (zero-based indexing)


Cut the bonds
-------------

   sep_mol.bonds.remove_bonds_between_atoms(iN, iN_C)
   sep_mol.bonds.remove_bonds_between_atoms(iO, iO_C)
   view(sep_mol, width=150, height=150, direction="tilt_pca3")

[image]


Increase the distance
---------------------

Now increase the distance. When you change the distance, all distances
within the separate molecules stay the same.

   sep_mol.set_distance(iN, iN_C, 2.5)  # angstrom
   view(sep_mol, width=150, height=150, direction="tilt_pca3")

[image]


Preoptimize the separated and ring systems
------------------------------------------

   from scm.plams import Settings, AMSJob, view

   preopt_s = Settings()
   preopt_s.input.ams.Task = "GeometryOptimization"
   preopt_s.input.DFTB.Model = "GFN1-xTB"
   preopt_ring_job = AMSJob(settings=preopt_s, molecule=ring_mol, name="preopt_ring_mol")
   preopt_ring_job.run()
   preopt_sep_job = AMSJob(settings=preopt_s, molecule=sep_mol, name="preopt_sep_mol")
   preopt_sep_job.run();

   [28.02|09:33:25] JOB preopt_ring_mol STARTED
   [28.02|09:33:25] JOB preopt_ring_mol RUNNING
   [28.02|09:33:26] JOB preopt_ring_mol FINISHED
   [28.02|09:33:26] JOB preopt_ring_mol SUCCESSFUL
   [28.02|09:33:26] JOB preopt_sep_mol STARTED
   [28.02|09:33:26] JOB preopt_sep_mol RUNNING
   [28.02|09:33:28] JOB preopt_sep_mol FINISHED
   [28.02|09:33:28] JOB preopt_sep_mol SUCCESSFUL

   <scm.plams.interfaces.adfsuite.ams.AMSResults at 0x71f14270a190>


Final molecule results
======================

   opt_ring_mol = preopt_ring_job.results.get_main_system()
   view(opt_ring_mol, width=150, height=150, direction="tilt_pca3")

[image]

   opt_sep_mol = preopt_sep_job.results.get_main_system()
   view(opt_sep_mol, width=150, height=150, direction="tilt_pca3")

[image]


Copy-pasteable system blocks
============================

When running calculations, bonds are really only used in ForceField
calculations. Otherwise, you do not need any bonding information (it
is not used). So it can be a bit misleading to store it or print it.
Therefore, we here remove all the bonds before printing.

   opt_ring_mol.bonds.clear_bonds()
   print(opt_ring_mol)

   System
      Atoms
         C 0.8176518265162523 0.006185364744877703 0.11652129166388352
         C -0.5531099061901967 -0.6363379504647356 -0.08529176454211394
         O -1.4164729952124286 0.48607450912407413 -0.19677333330677277
         N -0.6859030532374636 1.632537947036214 -0.08098563281192021
         N 0.5074709667559635 1.4283821326486141 0.08689807666792455
         H 1.2798721321383 -0.22791534308315498 1.0777663111036928
         H 1.5305100901334066 -0.2099050616627963 -0.6819652771198405
         H -0.866292408825807 -1.2485848049850283 0.7685478972017538
         H -0.6137266520780276 -1.2304367933580558 -1.0047175688566132
      End
   End

   opt_sep_mol.bonds.clear_bonds()
   print(opt_sep_mol)

   System
      Atoms
         C 1.0065881512963437 -1.417955004240563 0.12885420257373528
         C -0.08873018725926754 -2.1350000996605933 -0.03448649693154524
         O -1.5081108856471217 1.002080375323799 -0.2045570136880532
         N -0.4803377947492119 1.514182931120784 -0.052917508314934095
         N 0.4992559331764258 2.001851603909277 0.09160944741412266
         H 1.3565548327663348 -1.1143852514314225 1.1071131753303691
         H 1.6148659535066119 -1.095826278535031 -0.7065218137603257
         H -0.6985763255773102 -2.454234856430313 0.8007371044735704
         H -0.4402537210929026 -2.435675244924285 -1.0129788309893883
      End
   End


See also
========


Related examples
----------------

* Nudged Elastic Band (NEB) for gasphase reaction with DFTB

* Automated Transition-State Workflow for a Diels-Alder Addition

* Band structure with DFTB (SCC-DFTB, GFN1-xTB)


Related tutorials
-----------------

* NEB tutorial

* Building molecules tutorial


Related documentation
---------------------

* NEB task

* ChemicalSystem geometry manipulation


Python Script
=============

   #!/usr/bin/env python
   # coding: utf-8

   # ## Create the systems
   #
   # ### Create the ring molecule from SMILES
   #
   # First, we create ``ring_mol`` from the SMILES code, and then we cut the O-C and N-C bonds and increase the distance between the fragments. We need to first cut the bonds before increasing the distance, since ``ChemicalSystem.set_distance`` requires separate molecule fragments to work on.

   from scm.plams import view
   from scm.base import ChemicalSystem

   # from_smiles creates bonds between atoms
   ring_mol = ChemicalSystem.from_smiles("C1CON=N1")
   view(
       ring_mol,
       width=350,
       height=350,
       direction="tilt_pca3",
       show_atom_labels=True,
       atom_label_type="Name",
   )


   # ### Identify which bonds to cut
   #
   # It's easiest if you know the indices of the bonded atoms. Then you can just cut those bonds, e.g. ``ring_mol.bonds.remove_bonds_between_atoms(0, 4)``. That would remove the N-C bond in the above picture (there the indices start with 1, but when calling Python methods they start from 0)
   #
   # That is fine and the most straightforward way to cut the bonds. However, it is difficult to automate because the atomic indices may change for different structures.
   #
   # Here we show a way to programmatically identify "an N atom that is bonded to another N and C", and "an O atom that is bonded to an N and C". Then we cut those corresponding bonds.

   from typing import List, Tuple, Iterable, Optional


   def find_first_matching(
       cs: ChemicalSystem, elem: str, neighbors: Iterable[str]
   ) -> Optional[Tuple[int, List[int], List[str]]]:
       """Finds the first atom with element ``elem`` and neighbors ``neighbors``.

       Returns a 3-tuple: (index_of_atom, list_of_neighbor_indices, list_of_neighbor_elements)

       In the return value the list_of_neighbor_indices and list_of_neighbor_elements
       have the same length and correspond to each other.
       """

       sorted_target_neighbors = sorted(neighbors)
       for i, at in enumerate(cs):
           if at.symbol != elem:
               continue
           my_neighbors = list(cs.bonds.get_bonded_atoms(i))
           my_neighbors_elems = [cs.atoms[j].symbol for j in my_neighbors]
           if sorted(my_neighbors_elems) == sorted_target_neighbors:
               return i, my_neighbors, my_neighbors_elems

       return None


   sep_mol = ring_mol.copy()

   # N bonded to N and C
   iN, N_neighs, N_neighs_elems = find_first_matching(sep_mol, "N", ["N", "C"])
   iN_C = N_neighs[N_neighs_elems.index("C")]  # index of neighboring C atom

   # O bonded to N and C
   iO, O_neighs, O_neighs_elems = find_first_matching(sep_mol, "O", ["N", "C"])
   iO_C = O_neighs[O_neighs_elems.index("C")]  # index of neighboring C atom


   print(f"I will cut the bonds between atoms {iN}-{iN_C} and {iO}-{iO_C} (zero-based indexing)")


   # ### Cut the bonds

   sep_mol.bonds.remove_bonds_between_atoms(iN, iN_C)
   sep_mol.bonds.remove_bonds_between_atoms(iO, iO_C)
   view(sep_mol, width=150, height=150, direction="tilt_pca3", picture_path="picture1.png")


   # ### Increase the distance
   #
   # Now increase the distance. When you change the distance, all distances within the separate molecules stay the same.

   sep_mol.set_distance(iN, iN_C, 2.5)  # angstrom
   view(sep_mol, width=150, height=150, direction="tilt_pca3", picture_path="picture2.png")


   # ### Preoptimize the separated and ring systems

   from scm.plams import Settings, AMSJob, view

   preopt_s = Settings()
   preopt_s.input.ams.Task = "GeometryOptimization"
   preopt_s.input.DFTB.Model = "GFN1-xTB"
   preopt_ring_job = AMSJob(settings=preopt_s, molecule=ring_mol, name="preopt_ring_mol")
   preopt_ring_job.run()
   preopt_sep_job = AMSJob(settings=preopt_s, molecule=sep_mol, name="preopt_sep_mol")
   preopt_sep_job.run()
   # ## Final molecule results

   opt_ring_mol = preopt_ring_job.results.get_main_system()
   view(opt_ring_mol, width=150, height=150, direction="tilt_pca3", picture_path="picture3.png")


   opt_sep_mol = preopt_sep_job.results.get_main_system()
   view(opt_sep_mol, width=150, height=150, direction="tilt_pca3", picture_path="picture4.png")


   # ## Copy-pasteable system blocks
   #
   # When running calculations, bonds are really only used in ForceField calculations. Otherwise, you do not need any bonding information (it is not used). So it can be a bit misleading to store it or print it. Therefore, we here remove all the bonds before printing.

   opt_ring_mol.bonds.clear_bonds()
   print(opt_ring_mol)


   opt_sep_mol.bonds.clear_bonds()
   print(opt_sep_mol)
