How to make a python code that can read a .xyz file and find distance, angle and dihedral angle between atoms?

Question

I want to make a python script that will load an xyz file. From the xyz parameters, I need to find the distance between atoms, angle and dihedral between atoms.
If I have a xyz coordinate in this matrix form:

4
test
O 0.0 0.0 0.0
H 1.0 0.0 0.0
H 0.0 1.0 0.0

Then I want to find the bond parameters by simply using pythagoras formula using only the indices of the matrix. If I need to find the OH distance, I can give the indices and the value will come up.
Now this is easy for a small system, It will get tiresome for a large one. This is why I want to use a small script, that can load a xyz file, set index and allow user to find the bond parameters.

Answer 1

Assuming that you know how to calculate the dihedral angle algebraically, this python 3 code can be used. The data is from the xray structure of insulin 2INS. The angle is 136.7 degrees.

# Dihedral angles
phecd1 = np.array([ -7.497, 20.877, -0.865]) # xyz from 2INS.pdb  atom 349
phecd2 = np.array([ -5.644, 21.663,  0.493]) # atom 350 
tyrcd1 = np.array([ -4.106, 13.747,  5.282]) # atom 360  
tyrcd2 = np.array([ -1.869, 14.466,  4.723]) # atom 361  
vec_a = phecd2 - phecd1 
vec_b = tyrcd1 - phecd2  
vec_c = tyrcd2 - tyrcd1 

m = np.cross(vec_a,vec_b)
n = np.cross(vec_b,vec_c)
psi = np.sign(np.dot(vec_a,n))*np.arccos(np.dot(n,m)/( np.sqrt(np.dot(n,n)) *np.sqrt(np.dot(m,m)) ) )
print('{:s}{:6.2f}{:s}'.format('dihedral angle =',psi*180/np.pi,' degrees') )

Answer 2

It helps to know that in computational chemistry the collection of distance, angle and dihedral coordinates is collectively known as "z-matrix coordinates".

Unless you care to reinvent the wheel (which is not a bad exercise), it suffices to perform a web search for "python internal to zmatrix coordinates". For instance, the "chemcoord" python library claims to be capable of performing the desired transformation. Tutorials are available.

Answer 3

Here is an example for Open Babel 3 in Python 3.8. It's not really that pythonic, but you can make use of powerful features like SMARTS pattern and have the support for pretty much any chemical file format.

from openbabel import openbabel as ob

obconv = ob.OBConversion()
obconv.SetInFormat('xyz')
mol = ob.OBMol()
obconv.ReadFile(mol, 'test.xyz')
obsmarts = ob.OBSmartsPattern()
# Get O-H bonds
obsmarts.Init('O[H]')
if obsmarts.Match(mol):
    for match in obsmarts.GetUMapList():
        bond = mol.GetBond(*match)
        print(f'R{match} =', bond.GetLength())
# Get H-O-H angles
obsmarts.Init('[H]O[H]')
if obsmarts.Match(mol):
    for match in obsmarts.GetUMapList():
        print(f'A{match} =', mol.GetAngle(*(mol.GetAtom(i) for i in match)))

Answer 4

GeomeTRIC and pyGSM have this capability for internal coordinates

In pyGSM this looks like

from pygsm.utilities.manage_xyz import read_xyz,xyz_to_np
from pygsm.coordinate_systems.slots import Distance,Angle,Dihedral

bond = Distance(0,1)
angle = Angle(0,1,2)
dihedral=Dihedral(0,1,2,3)

geom = read_xyz('ethylene.xyz')
xyz = xyz_to_np(geom)

bond_value = bond.value(xyz)
angle_value = angle.value(xyz)
dihedral_value = dihedral.value(xyz)

print(bond_value,angle_value,dihedral_value)