Generate structure using atomic coordinates

Question

I am trying to write a Python/C++ code that reads in a *.xyz file and varies the interatomic distances, angles etc. between distinct atoms in a molecule. In other words, a complete new structure with varied bond lengths should be generated. It's clear to me how to calculate the distance in cartesian coordinates between two points (atoms) P and Q.

For a better understanding I used this Methanol structure provided by the Avogadro suite:

6
H          0.84754        0.03474        1.03453
C          0.35048        0.00675        0.06084
H          0.63507        0.89276       -0.52007
H          0.66294       -0.89338       -0.48283
O         -1.01083       -0.00823        0.36439
H         -1.48520       -0.03263       -0.45685

However, when I manually vary the bond length between the atoms C-O by 0.1 Angstrom in Avogadro, e.g., in Methanol I have difficulties understanding how to find the new cartesian coordinates that are translated/transformed in the newly generated structure:

6
H          0.94513        0.03581        1.01277
C          0.44807        0.00782        0.03908
H          0.73266        0.89383       -0.54183
H          0.76053       -0.89231       -0.50459
O         -1.01083       -0.00823        0.36439
H         -1.48520       -0.03263       -0.45685

Could someone explain to me the mathematical step for generating this new set of cartesian coordinates or direct me to an example?

By looking into the Header files of Avogadro, I can tell that computationally one needs to or can use the API, i.e., the main module of Eigen which allows dense matrix and array manipulation.

I am grateful for any help!

Answer 1

To understand the new coordinates after altering the bond length of the C-O bonds you can use the below approach:

Let us first visualize the positions of all the atoms before and after altering the bond length:

The atoms with * are the atoms at the new coordinates after bond altering.

So imagine a vector from the O atom to the C atom. We would like to get the new coordinates of the C atom after increasing the distance between the current C atom and the O atom by 0.1 Angstrom.

This can be achieved by the following formula:

D = dist_O_and_C_before_bond_alter + 0.1
C_atom_new_pos_vector = O_vector_atom_before_alter + D * OC_unit_vector

So to calculate C_atom_new_pos_vector we need to calculate three other elements: {2 scalars: (dist_O_and_C_before_bond_alter, D) and 1 vector: OC_unit_vector}

The code to calculate all these elements is shown below in python:

atom_coords_before_bond_alter = np.array([
[0.84754,0.03474,1.03453], # H
[0.35048,0.00675,0.06084], # C
[0.63507,0.89276,-0.52007], # H
[0.66294,-0.89338,-0.48283], # H
[-1.01083,-0.00823,0.36439], # O
[-1.48520,-0.03263,-0.45685] #H
])

def get_dist(atom_1_coords, atom_2_coords):
  x_1 = atom_1_coords[0]
  y_1 = atom_1_coords[1]
  z_1 = atom_1_coords[2]

  x_2 = atom_2_coords[0]
  y_2 = atom_2_coords[1]
  z_2 = atom_2_coords[2]

  d = ( (x_1-x_2)**2 + (y_1-y_2)**2+ (z_1-z_2)**2 )**0.5
  return d

dist_O_and_C_before_bond_alter = get_dist(atom_coords_before_bond_alter[1],atom_coords_before_bond_alter[4])

C_vector_atom_before_alter = atom_coords_before_bond_alter[1]
O_vector_atom_before_alter = atom_coords_before_bond_alter[4]
OC_dir_vector = C_vector_atom_before_alter - O_vector_atom_before_alter

OC_vector_magnitude = np.sqrt(np.sum([comp**2 for comp in OC_dir_vector]))
OC_unit_vector = OC_dir_vector/OC_vector_magnitude
D = dist_O_and_C_before_bond_alter + 0.1

C_atom_new_pos_vector = O_vector_atom_before_alter + D * OC_unit_vector

Let us now print C_atom_new_pos_vector; This outputs >> array([0.44807731, 0.00782397, 0.03907739]) which matches the vector that you get from Avagadro for the new position of the C atom after bond altering

So similarly, you can calculate the new positions of the other atoms to understand the new coordinates of the other atoms.

Answer 2

The easiest way to do this IMHO is through Open Babel / Pybel, particularly if you want to set a lot of bonds.

BTW, reading from XYZ, the bond orders are interpreted, but not always exactly.

In Avogadro v1, this is basically the code for setting bond lengths. In Avogadro v2, it doesn't use Open Babel directly, but the code is quite similar.

import sys
import os

from openbabel import pybel
# don't need to reimport openbabel
ob = pybel.ob

# syntax:
# set-bondlength.py file.xyz 

# repeat through all the files on the command-line
# we can change this to use the glob module as well
#  e.g., find all the files in a set of folders
for argument in sys.argv[1:]:
    filename, extension = os.path.splitext(argument)

    for i in range(len(filename), 0, -1):
        if filename[i-1].isalpha():
            field = filename[i:]
            break

    # read all the molecules from the supplied file
    for mol in pybel.readfile(extension[1:], argument):

        for bond in ob.OBMolBondIter(mol.OBMol):
            beginAtom = bond.GetBeginAtom()
            endAtom = bond.GetEndAtom()

            # check if the bond is a C-C single bond
            if beginAtom.GetAtomicNum() == 6 \
                    and endAtom.GetAtomicNum() == 6 \
                    and bond.GetBondOrder() == 1:

                # set the bond length to 1.54 Angstrom
                bond.SetLength(1.54)
        
        # write the molecule to a new file
        mol.write("xyz", filename + "_%s.xyz" % field, overwrite=True)