9
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ You mean you want us to write a script? Do you mind if I ask: Can you program in Python? Have you tried yourself? Have you tried to Google if such library/module exist? I think this is a HW question, with no effort. $\endgroup$
    – Greg
    Mar 13 '21 at 5:38
  • $\begingroup$ I am sorry I this question felt like this. I am new to python. I found examples where I could use the xyz coordinates as tuples but I had to write them explicitly in the code. I am still doing that. However, I am also trying if the code can read the coordinates and treat the coordinates as Xij, Yij and Zij element so then I can take any .xyz file and use the code to find the distance. Again I am sorry if my question felt like this. $\endgroup$
    – Prasanta Bandyopadhyay
    Mar 13 '21 at 11:30
13
$\begingroup$

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') )

dihedral angle

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you @porphyrin. One small favour, did you import this line manually from the .xyz file? phecd1 = np.array([ -7.497, 20.877, -0.865]) If no, how can one do that? $\endgroup$
    – Prasanta Bandyopadhyay
    Mar 5 '21 at 13:59
  • 2
    $\begingroup$ This was just an example. The easiest way if you are using protein data bank .pdb files is to read in the file (see examples in python /numpy info, eg readlines, genfromtxt) as text then look for ATOM or HETATOM as the first word then select x,y,z from that row and put into an array. If its other file types you many need to ignore first so many lines and again select coordinates and type of atom if that is what you need. $\endgroup$
    – porphyrin
    Mar 5 '21 at 15:27
  • $\begingroup$ Suppose you want to drag one atom along the z-axis by 1 A untill the vec_a becomes 10 A. How can you do that? I guess the way should be using the while loop. But I can not use the while loop for a certain element of the array. $\endgroup$
    – Prasanta Bandyopadhyay
    Mar 10 '21 at 4:34
  • 1
    $\begingroup$ If you want to do this to say Cb on each residue then you will make a loop over all atoms and inside this loop look for atoms of type Cb and make the change. $\endgroup$
    – porphyrin
    Mar 10 '21 at 7:45
12
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank You Buck Thorn. I am looking at the tutorial and chemcoord.readthedocs.io/en/v2.0.4. This will be helpful. It was just, that I wanted to sort the matrix as ijk th element and perform further calculations. $\endgroup$
    – Prasanta Bandyopadhyay
    Mar 3 '21 at 11:03
  • $\begingroup$ There are ofc many more libraries to pick from. Cheminformatics toolkits that offer Python wrappers like Open Babel or RDKit would also get the job done and offer a wide range of other functionalities in addition. $\endgroup$
    – Kexanone
    Mar 13 '21 at 2:06
8
$\begingroup$

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)))
Improve this answer
$\endgroup$
6
$\begingroup$

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)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.