Breaking bonds in a ring molecule and modifying distances for NEB

Building structures Technical DFTB NEB

Downloads: Notebook | Script ?

Requires: AMS2026 or later

Related examples
Related tutorials
Related documentation
../_images/building_structures_cut_ring_molecule_for_neb_14_0_0668afe5.png

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",
)
../_images/building_structures_cut_ring_molecule_for_neb_1_0_cb276282.png

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")
../_images/building_structures_cut_ring_molecule_for_neb_7_0_cea480b0.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")
../_images/building_structures_cut_ring_molecule_for_neb_9_0_08146a85.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();

[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")
../_images/building_structures_cut_ring_molecule_for_neb_13_0_7ee16568.png 
opt_sep_mol = preopt_sep_job.results.get_main_system()
view(opt_sep_mol, width=150, height=150, direction="tilt_pca3")
../_images/building_structures_cut_ring_molecule_for_neb_14_0_0668afe5.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)

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

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)